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NOMENCLATURE 


a  Dimensionless  wavelength 

AR  Aspect  ratio 

B  Parameter  in  equation  of  state  for  vapor 
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W  Cell  dimensions 
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THERMAL  CONVECTION  IN  SNOW 


D.J.  Powers,  S.C.  Colbeck  and  K.  O’Neill 

INTRODUCTION 

The  size  and  shape  of  grains  and  the  density  of  snow  greatly  affect  all  of  its  material  properties. 
Snow  covers  of  similar  density  and  grain  size  may  have  widely  different  crystal  shapes,  ranging  from 
rounded  to  sharply  angular.  The  range  over  which  the  physical  properties  vary  is  quite  large. 

It  is  widely  known  that  when  a  low  density  dry  snow  is  subjected  to  a  temperature  gradient  of 
at  least  0.1°C/cm,  many  of  the  grains  will  develop  facets.  A  dramatic  loss  of  strength  usually  ac¬ 
companies  this  rapid  recrystallization  and  this  strength  loss  is  commonly  thought  to  be  a  factor  in 
the  release  of  avalanches.  Efforts  to  model  crystal  growth  in  snow  have  so  far  focused  solely  on 
vapor  diffusion  among  the  crystals.  However,  the  temperature  gradients  that  drive  the  diffusion 
of  water  vapor  also  establish  air  density  gradients  that  may  lead  to  convective  flows  of  air.  These 
flows,  if  they  occur,  would  certainly  have  a  large  impact  on  the  flux  of  heat  and  vapor  and  there¬ 
fore  might  affect  the  metamorphism  of  snow. 

Snow  metamorphism 

Dry  snow  metamorphism  is  often  classified  as  equitemperature  or  temperature  gradient.  The 
first  term  refers  to  conditions  under  which  grains  become  more  rounded  and  sintering  processes 
increase  the  size  of  bonds  between  grains.  It  is  now  recognized  that  these  processes  do  not  ever 
occur  under  isothermal  conditions,  as  the  name  would  imply.  Furthermore,  when  nature  imposes 
even  a  slight  gradient  on  the  snow  cover  the  processes  are  greatly  accelerated.  Colbeck  (1982)  re¬ 
fers  to  the  rounded  grains  as  the  equilibrium  form,  which  is  probably  a  better  description  than 
equitemperature. 

Temperature  gradient  metamorphism  is  characterized  by  the  growth  of  large  angular  or  faceted 
crystals.  Actually,  temperature  gradient  metamorphism  is  also  a  misnomer,  since  large  faceted 
crystals  occur  only  when  large  temperature  gradients  are  applied  to  low  density  snow.  A  loss  of 
strength  frequently  accompanies  the  growth  of  large  faceted  grains,  especially  when  hollow  depth 
hoar  crystals  form.  These  grow  most  readily  in  low  density  snows  subjected  to  very  large  temper¬ 
ature  gradients.  More  solid  faceted  crystals  grow  at  lower  temperature  gradients,  or  in  higher  den¬ 
sity  snows  (Akitaya  1974,  Marbouty  1980). 

Mass  transfer  by  diffusion  in  snow 

Here  we  are  primarily  concerned  with  metamorphism  driven  by  strong  temperature  gradients 
because  the  temperature  gradients  that  drive  depth  hoar  growth  could  also  drive  the  convective 
flows  that  are  the  subject  of  the  present  work. 

Metamorphism  in  dry  snow  occurs  by  the  transport  of  water  vapor.  Vapor  diffuses  from  areas 
of  high  vapor  density  to  low  vapor  density,  or  it  may  be  carried  by  convective  currents  of  air.  In 


this  section,  we  discuss  only  the  diffusion  process.  Vapor  pressure  gradients  exist  in  snow  because 
of  temperature  gradients  or  curvature  differences,  or  both.  Vapor  pressure  is  higher  over  warmer 
ice  and  over  convexly  curved  ice.  When  temperature  gradients  are  low,  curvature  effects  are  im¬ 
portant,  and  concave  sections  will  begin  to  fill  in  by  vapor  deposition.  However,  at  the  large  val¬ 
ues  of  temperature  gradient  where  depth  hoar  grows,  curvature  effects  are  negligible,  and  the 
transport  of  vapor  is  governed  by  the  temperature  field.  Strong  gradients  occur  for  prolonged 
periods  most  commonly  in  cold  climates  with  shallow  snow  covers  (less  than  100  cm). 

If  we  think  of  a  vertical  coordinate  as  positive  upwards,  the  gradient  will  be  negative.  Hence¬ 
forth,  a  large  gradient  refers  to  a  large  absolute  value  of  the  negative  gradient. 

Early  attempts  to  model  depth  hoar  growth  (i.e„  Giddings  and  LaChapelle  1962)  used  a  one¬ 
dimensional  diffusion  equation  to  describe  the  flux  of  vapor  through  dry  snow: 


(1) 


where  D  is  a  diffusion  coefficient.  Giddings  and  LaChapelle  tried  to  calculate  crystal  growth  rates 
assuming  that  all  vapor  stopped  at  one  level.  This  assumption  is  generally  meaningless,  although 
eq  1  can  be  used  to  represent  mass  fluxes  in  one  dimension.  Experimental  data  from  Yosida  (1955) 
and  Yen  (1963)  show  a  rate  of  vapor  transport  four  or  five  times  greater  than  the  Giddings  and  La¬ 
Chapelle  model  predicts.  Yosida  formulated  his  results  in  an  equation  similar  to  eq  1,  using  an  ef¬ 
fective  diffusivity  coefficient  to  describe  the  flux  of  vapor.  He  found  De ^  equal  to  0.85  cmJ/s 
and  invariant  with  snow  density  over  the  range  of  0.08  <  ps  <  0.51  g/cm3  .  Yen  confirmed  these 
results,  although  he  tested  only  densities  in  the  range  of  0.38  <  ps<  0.48  g/cm3 .  Trabant  and 
Benson  (1972)  calculated  flux  rates  in  the  field  and  found  them  to  be  about  an  order  of  magni¬ 
tude  greater  than  they  calculated  using  eq  1 .  Although  they  did  not  calculate  an  effective  diffu¬ 
sivity,  their  results  indicate  an  effective  diffusivity  about  a  factor  of  2  greater  than  that  of  Yosida. 

The  reason  for  the  high  observed  flux  rates  is  not  entirely  clear.  Trabant  and  Benson  suggest 
that  thermal  convection  may  be  important,  and  this  may  certainly  be  part  of  the  answer.  In  Yo- 
sida’s  experiments  the  geometry  would  make  convection  unlikely.  The  common  observation  of 
rapidly  growing  depth  hoar  in  shallow  layers  of  snow  where  convection  is  unlikely  also  points  to 
the  need  for  some  other  explanation  for  the  high  flux  rates. 

Yosida  supposed  that  vapor  transport  occurred  by  a  “hand-to-hand”  mechanism,  i.e.,  that  mass 
was  transferred  by  diffusion  from  grain  to  grain  and  not  along  vertical  air  channels.  Since  ice  is  a 
much  better  conductor  than  air,  temperature  gradients  in  the  air  spaces  may  be  much  higher  than 
the  average  gradient  in  the  snow,  and  thus  local  fluxes  may  be  much  higher  than  that  predicted  on 
the  basis  of  an  average  temperature  gradient.  This  is  part  of  the  explanation  for  the  high  values  of 
Delj.  Colbeck  ( 1983)  quantified  this  effect  in  an  effort  to  model  grain  growth.  His  model  did  show 
reasonable  growth  rates,  although  he  was  forced  to  arbitrary  assumptions  about  snow  stereology  be¬ 
cause  of  a  lack  of  good  data.  His  model  has  not  been  applied  to  calculate  macroscopic  mass  flux 
rates,  and  thus  it  is  not  yet  certain  just  what  the  quantitative  effects  of  hand-to-hand  transfer  are 
on  those  mass  fluxes. 

Colbeck’s  model  showed  that  grain  growth  is  primarily  attributable  to  a  coupling  between  verti¬ 
cally  aligned  particles,  just  as  Yosida  suggested.  The  rate  of  mass  movement  between  the  particles 
is  equal  to  the  rate  of  mass  gain  by  the  growing  particle.  Thus  grains  do  not  grow  because  of  a 
macroscopic  redistribution  of  mass  between  regions  of  the  snow  cover  as  described  in  some  flux 
divergence  models.  The  particle-to-particle  mechanism  explains  the  experimental  observation  of 
rapid  grain  growth  without  significant  snow  density  change  (Marbouty  1980).  The  relatively  in¬ 
tense  vapor  transport  between  coupled  grains  is  included  in  the  overall  macroscopic  vapor  flux  de¬ 
scribed  by  eq  1 ,  with  an  effective  diffusion  coefficient.  We  study  the  macroscopic  fluxes  and  their 
possible  enhancement  by  convection  in  order  to  see  how  the  convective  fluxes  affect  the  grain-scale 
vapor  transport  that  controls  snow  metamorphism. 


Heat  transfer 

Typically,  heat  transfer  through  snow  is  modeled  by  a  conduction  equation  of  the  form 


,  3  T 
Qz=~kHz 


(2) 


where  k  is  a  thermal  conductivity.  It  is  customary  to  use  an  effective  thermal  conductivity  to  ac¬ 
count  for  the  heat  conduction  through  both  the  fluid  and  solid  in  porous  media.  In  snow  we  must 
also  consider  the  latent  heat  carried  by  the  flux  of  vapor.  The  transport  of  heat  ascribable  to  vapor 
diffusion  is  given  by 


dp^dT 

qV2  L  Deff  gy, 


(3) 


where  L  is  the  latent  heat  of  sublimation. 

Since  this  flux  is  also  proportional  to  the  temperature  gradient,  all  contributions  to  heat  flow 
may  be  lumped  into  an  effective  thermal  conductivity,  as  described  by  eq  2.  Data  on  these  effec¬ 
tive  thermal  conductivities,  summarized  by  Mellor  (1977),  show  a  great  deal  of  scatter  for  any 
given  density.  In  part  this  may  be  attributable  to  the  temperature  dependence  of  the  vapor  term. 
The  experimental  data  of  Pitman  and  Zuckerman  (1967)  show  the  expected  trends  with  tempera¬ 
ture;  however,  the  magnitude  of  the  change  in  keff  with  temperature  that  they  observed  is  too 
great  to  be  explained  by  the  temperature  dependence  of  the  vapor  diffusion  term. 

A  basic  understanding  of  the  relative  contributions  of  conduction  through  ice,  of  conduction 
through  air  and  of  vapor  diffusion  to  the  total  transport  of  heat  is  lacking.  Generally,  it  is  assumed 
that  conduction  through  the  air  is  negligible,  since  the  thermal  conductivity  of  air  is  100  times  less 
than  that  of  ice.  However,  the  aforementioned  notion  of  exaggerated  gradients  in  the  voids  be¬ 
tween  particles,  coupled  with  the  high  porosity  of  snow,  leads  us  to  estimate  that  the  conduction 
of  heat  through  the  air  spaces  may  be  on  the  order  of  10%  of  that  through  the  ice  lattice. 

The  contribution  of  vapor  transport  can  be  estimated  from  eq  3,  defining  kv  as  the  vapor  diffu¬ 
sion  equivalent  of  thermal  conductivity,  or 


bpv 

kv  —  L  f^eff  0  jr  ' 


(4) 


where  kv  is  a  function  of  temperature,  as  illus¬ 
trated  in  Figure  1 .  At  0°C,  using  Yosida’s  De„ 
of  0.85  cmJ/s,  we  find  kv  =  2.4x  HT4  cal/cm 
s  °C.  Since  total  thermal  conductivities  of  less 
than  this  are  observed,  we  must  be  somewhat 
skeptical  of  either  Yosida’s  notion  that  Deff  is 
not  a  function  of  density,  or  of  the  thermal 
conductivity  data  (or  of  both).  It  is  reasonable 
that  Deff  is  a  function  of  density,  since  crystal 
growth  rates  are  strongly  dependent  on  snow 
density. 

In  summary,  for  lower  density  snows,  con¬ 
duction  through  air  is  probably  not  negligible, 
while  the  transfer  of  latent  heat  is  probably  very 
important.  In  higher  density  snows,  where  the 
total  thermal  conductivity  is  greater,  the  po- 
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Figure  1.  Contributions  of  vapor  diffusion 
to  heat  transport  in  terms  of  thermal  con¬ 
ductivity  kv  for  various  values  of  effective 
diffusivity  D. 


rosity  is  less  and  temperature  differences  among  neighboring  grains  are  less,  conduction  through 
the  ice  lattice  may  dominate  the  heat  transfer  process. 


BACKGROUND-POROUS  MEDIA 

This  section  examines  the  developments  that  have  taken  place  in  the  study  of  thermal  convec¬ 
tion  in  porous  media.  The  impetus  for  much  of  this  work  has  been  an  interest  in  petroleum  and 
geothermal  reservoirs.  However,  much  of  the  work  is  of  a  basic  nature,  and  may  be  applied  to 
snow.  The  aim  here  is  to  outline  the  fundamental  principles  that  govern  thermal  convective  flows 
in  porous  media  and  to  summarize  the  pertinent  results. 

Structure  of  thermal  convection 

Convective  tlows  occur  when  buoyancy  forces,  driven  by  density  gradients,  are  sufficient  to 
overcome  viscous  drag.  As  warm  air  rises,  cold  air  from  above  must  descend  for  continuity.  Ther¬ 
mal  convection  can  thus  be  thought  of  as  nature’s  way  of  trying  to  relieve  unstable  density  gradi¬ 
ents. 

Circulation  patterns  in  a  horizontal  layer  heated  from  below  appear  commonly  in  two  forms. 
One  form,  commonly  called  two-dimensional  convection,  consists  of  rolls  (see  Fig.  2)  that  are 
very  long  in  one  direction,  and  are  roughly  equal  in  the  other  two  dimensions.  The  other  form  is 
three-dimensional  convection,  in  which  the  cell  takes  a  hexagonal  pattern  when  viewed  from 
above.  The  cell  size  at  the  onset  of  convection  was  calculated  for  both  forms  by  Combarnous  and 
Bories  (1975)  (Fig.  2). 

When  the  porous  layer  is  confined  laterally,  the  dimensions  of  the  convective  cells  are  governed 
by  the  dimensions  of  the  container.  Beck  (1972),  and  Tewari  and  Torrance  (1981)  have  predict¬ 
ed  the  preferred  cell  size  at  the  onset  of  convection  as  a  function  of  container  dimensions  for  the 
closed-top  and  open-top  cases  respectively.  Cheng  (1978),  in  his  review  paper,  draws  two  general 
conclusions  about  the  results  of  Beck,  and  these  seem  to  apply  equally  well  to  the  results  of  Tew¬ 
ari.  First,  when  one  of  the  lateral  dimensions  is  less  than  the  vertical  dimension,  two  dimensional 
convection  is  the  rule.  Secondly,  the  number  and  direction  of  rolls  is  governed  by  the  tendency 
for  the  width  of  the  roll  to  be  nearly  equal  to  its  height.  We  can  also  say  that  when  the  upper 
boundary  is  permeable  the  lateral  dimensions  of  the  convective  cell  are  about  one  and  one-half 
times  greater  than  for  the  closed  top  case,  all  other  things  being  equal. 
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The  above  results  strictly  apply  only  at  the  onset  of  convection.  When  convection  is  “mild,” 
these  results  are  probably  a  good  approximation.  However,  as  convection  intensifies,  the  pre¬ 
ferred  cell  size  can  be  greatly  reduced.  Numerical  experiments  by  Straus  and  Schubert  (1979) 
show  that  for  random  initial  conditions,  either  two-  or  three-dimensional  convection  can  occur. 
The  mode  that  would  maximize  heat  transfer  is  not  necessarily  the  one  that  results. 

Rayleigh  number 

The  flow  of  fluids  through  a  porous  medium  is  described  by  Darcy’s  law,  which  is  an  empiri¬ 
cal  relation  developed  by  Darcy  from  his  classic  set  of  experiments.  Two  assumptions  are  made 
here:  that  the  fluid  is  Boussinesq,  meaning  that  variations  in  fluid  properties  are  negligible  except 
where  buoyancy  terms  are  involved,  and  that  it  is  incompressible,  which  implies  that  density  will 
not  vary  directly  as  a  result  of  the  flow. 

From  a  nondimensional  form  of  Darcy’s  law  (derived  more  completely  in  the  Modeling  section) 
comes  an  important  nondimensional  parameter  known  as  the  Rayleigh  number,  defined  as 
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Ra  =  — -  (5) 
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where  k  is  kmKpCp)f  and  p0  is  the  fluid  density  at  a  reference  temperature.  The  Rayleigh  num¬ 
ber,  along  with  the  boundary  conditions,  governs  both  the  onset  of  convection  and  its  intensity. 

The  definition  of  the  temperature  difference  in  the  Rayleigh  number  depends  on  the  thermal 
boundary  conditions  at  the  upper  and  lower  surfaces.  When  the  boundary  conditions  are  con¬ 
stant  temperatures  at  both  upper  and  lower  surfaces,  then  the  AT  is  just  the  difference  between 
those  temperatures.  However,  we  are  also  interested  in  lower  boundaries  with  a  constant  heat 
flux  (<7 ).  When  this  is  the  case,  we  define  the  temperature  difference  as  qH/km ,  which  is  the  tem¬ 
perature  difference  that  would  be  observed  if  convection  were  absent. 

Onset  problem 

Convection  takes  place  when  buoyancy  forces  are  sufficient  to  overcome  viscous  resistance. 

For  a  porous  medium,  this  happens  when  the  Rayleigh  number  is  sufficiently  large.  Die  critical 
value  of  the  Rayleigh  number  at  which  convection  begins  is  a  function  primarily  of  the  boundary 
conditions.  The  critical  Rayleigh  number  was  originally  calculated  for  a  horizontal  layer  with  iso¬ 
thermal.  impermeable  boundaries  by  Lapwood  (1948).  Nield  (1968)  has  applied  the  analy  sis  of 
Lapwood  to  other  boundary  conditions,  including  permeable  or  constant  heat  flux  boundaries. 
Some  of  his  results,  pertinent  to  the  present  problem,  are  listed  in  Table  1 ,  along  with  the  results 
of  the  original  problem  considered  by  Lapwood. 


Table  1.  Summary  for  the  onset  of  convection.  After 
Nield  (1968). 


Case 

Top  boundary 

Bottom  boundary  Racr 

acr 

1 

Constant  temperature. 

Constant  temperature,  39.5 

3.14 

impermeable 

impermeable 

2 

Constant  temperature, 

Constant  temperature.  27.1 

2.3  3 

permeable 

impermeable 

3 

Constant  temperature. 

Constant  flux,  27.1 

2.33 

impermeable 

impermeable 

4 

Constant  temperature. 

Constant  flux.  17.7 

1.75 

permeable 

impermeable 

uses  eq  33  for  the  energy  equation,  and  thus  accounts  for  latent  heat,  while  the  other  uses  eq  3 1 , 
and  does  not  account  for  latent  heat  effects.  The  codes  are  designated  as  CONVAP  and  CONVEC 
respectively.  The  programs  differ  only  in  the  energy  equation  used,  and  thus  Figure  9  applies  to 
both  codes. 

Tire  initial  temperature  profile  varied  linearly  in  z  and  was  constant  in  the  x  direction.  The 
stream  function  was  set  equal  to  zero  over  the  entire  grid,  corresponding  to  a  no-flow  situation. 

A  temperature  perturbation  was  introduced  prior  to  the  first  stream  function  calculation  to  start 
the  Hows.  The  location  and  magnitude  of  this  perturbation  affects  the  transient  response  but  not 
the  final  steady  state. 

Typically,  most  of  the  computing  time  was  used  to  obtain  the  iterative  solution  of  the  equa¬ 
tion  of  motion.  This  equation  was  therefore  overrelaxed  (Pinder  and  Gray  1977)  to  speed  con¬ 
vergence.  In  essence,  what  the  overrelaxation  does  is  to  weight  the  value  of  t/zfj1  to  make  the 
equation  more  or  less  implicit.  It  was  found  that  an  optimum  value  of  the  relaxation  coefficient 
was  1 .7,  where  2.0  is  fully  implicit  and  0.0  is  fully  explicit.  The  convergence  of  the  solution  at 
a  time  step  was  adequate  when  for  all  i  and  / 
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where  &max  and  4>min  are  the  maximum  and  minimum  values  of  the  stream  function  calculated 
during  iteration. 

Several  criteria  were  used  to  determine  if  the  solution  was  at  a  steady  state.  The  most  reliable 
was  to  monitor  energy  conservation.  At  steady  state,  all  heat  flowing  in  must  flow  out.  Since 
the  lateral  boundaries  are  adiabatic,  this  implies  that  the  heat  flux  at  the  bottom  be  equivalent  to 
the  heat  flux  at  the  top.  In  nondimensional  form,  the  heat  flux  is  the  Nusselt  number.  The  con¬ 
servation  of  energy  was  monitored  by  observing  the  ratio  |(Nurop  -  Nuftof)/Nuhof|.  The  equa¬ 
tions  used  for  calculating  the  Nusselt  numbers  at  the  boundaries  are  derived  in  Appendix  A. 

For  each  run,  the  energy  balance  ratio  would  reach  a  minimum  value  because  of  accuracy  lim¬ 
itations  of  the  solution.  Better  conservation  (a  lower  minimum  balance)  could  be  achieved  by  re¬ 
fining  the  grid.  Choice  of  the  time  step  did  not  affect  the  final  steady  state  or  the  energy  conser¬ 
vation  achieved.  Thus  the  final  steady  state  was  independent  of  the  path  followed  to  obtain  it. 

Choice  of  the  proper  time  step  was  important  in  the  speed  with  which  the  solution  converged. 
Too  small  a  time  step  resulted  in  smooth  but  very  slow  convergence.  Too  large  a  time  step  would 
produce  an  oscillating  solution  that  would  overshoot  the  steady  state,  thus  also  requiring  a  long 
solution  time.  If  the  time  step  was  much  too  large,  the  solution  would  not  converge  to  a  steady 
state.  From  experience  we  found  that  the  optimum  time  step  size  was  usually  between  0.01  and 
0.05.  Grid  spacing  size  did  not  seem  to  affect  the  choice  of  Ar,  but  the  value  of  the  Rayleigh 
number  did.  It  was  found  to  be  necessary  to  use  a  smaller  time  step  as  the  Rayleigh  number  in¬ 
creased. 

Verification  of  the  model 

To  test  the  model,  runs  were  made  that  could  be  compared  directly  to  results  in  the  literature. 
Most  of  the  results  in  the  literature  are  for  the  basic  case  of  convection  in  a  horizontal  layer  bound¬ 
ed  by  isothermal,  impermeable  surfaces.  Figure  3  shows  some  of  the  results  for  this  case.  Experi¬ 
mental  results  (compiled  by  Combarnous  and  Bories)  are  shown,  along  with  analytical  results  ot 
Palm  et  al.  ( 1972)  and  Combarnous  and  Bories  (1975).  The  simple  relation  deduced  by  Elder 
(1967). 

Nu  =  Ra/Raff.  (43) 
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hxh, 

accumulation  =  iT^F'  -  T\  j)— — 


conduction 


ants 


convection  =  -  (uT)j+K  h,At  +  [(wT),.^  ~(wT)f+v J  —  At 
Rearranging  and  combining,  we  find  the  final  form  to  be 

(T*+'-T<\  JT2j-Tu\  TlJ+i-2TIJ+Tl<hl 

\  Ar  a /‘hn 

(42) 

2(uT)l  f ,/;  RwD^vj,  -(w7>.yi'[ 

K  "L  hz  J- 

While  the  correct  form  of  the  conduction  terms  could  be  found  front  a  shadow  node  technique 
similar  to  the  one  employed  earlier,  the  correct  form  of  the  convection  term  is  much  more  diffi¬ 
cult  to  derive. 

Numerical  solution 

An  Alternating  Direction  Implicit  (ADI)  scheme  is  used  to  evaluate  the  algebraic  equations 
(Pinder  and  Gray  1977).  The  method  sweeps  the 
grid  line  by  line,  first  in  the  x  direction,  then  in  the 

z  direction.  The  scheme  is  implicit  in  that  the  equa-  f  ~  ] 

.  Read  constants  I 

tions  written  about  each  node  are  solved  simultane-  _ , _ j 

ously  for  each  successive  line.  Thus  each  algebraic- 
equation  contains  more  than  one  unknown.  When 
solving  the  equations  for  a  line  swept  in  the*  direc¬ 
tion,  the  derivatives  in  *  are  evaluated  at  the  time 
step  level  d+1  and  derivatives  in  z  are  evaluated  at 
the  time  step  level  ?.  The  alternatives  to  a  line- 
implicit  scheme  are  to  solve  the  equations  explicitly, 
where  new  values  are  calculated  on  a  node  by  node 
basis,  or  to  solve  the  equations  directly,  in  which 
case  the  values  for  the  entire  grid  are  calculated  si¬ 
multaneously.  The  ADI  method  is  more  stable  than 
the  explicit  method,  thus  allowing  much  larger  time 
steps  to  be  used.  The  advantage  over  direct  methods 
is  that  the  matrix  of  coefficients  is  not  zero  only  on 
three  main  diagonals,  and  thus  the  system  of  equa¬ 
tions  may  be  solved  simply  and  efficiently  without 
evaluating  all  of  the  zeros.  An  explicit  form  of  this 
program  was  originally  developed,  and  was  found  to 
be  much  slower  and  more  costly. 

The  computing  flow  chart  is  given  in  F  igure  9,  Figure  9.  Flow  diagram  for  computer 

Two  FORTRAN  programs  were  developed:  One  codes  CONVAP and  CONVF.C. 


u 


18 


ary.  Beier  et  al.  (1983)  account  for  variations  along  the  boundary  and  find  that  the  second  order 
terms  are  of  0(/i3 )  and  0(/i4 )  for  the  convective  and  conductive  terms  respectively.  Thus  each  is 
less  than  the  first  order  truncation  error  term  by  0 (h2)  and  can  be  neglected  for  reasonably  fine 
grids. 

The  complete  energy  equation  is  then 


Ar  "  h2x  h2z 


(40) 


The  terms  on  the  right-hand  side  may  be  evaluated  at  either  the  £  or  £+1  time  step  level,  depend¬ 
ing  on  the  solution  scheme  used.  This  same  form  can  be  obtained  by  applying  central  differences 
to  the  conduction  terms,  first  order  upwind  differences  to  the  convection  term,  and  a  first  order 
forward  difference  to  the  transient  term. 

When  the  effects  of  phase  change  are  considered  (eq  32),  the  finite  differences  used  follow  the 
same  basic  concepts.  Convection  of  vapor  is  treated  in  a  similar  way  to  the  convection  of  air. 
Diffusion  of  vapor  is  treated  in  a  manner  similar  to  that  of  conduction  of  heat.  The  discrete  form 
of  eq  3 1  (again  with  a  false  transient)  is 


AWWr,)-*wr,-rM) 

Ar  ‘  hi 

|^(^i-r/)-^2/.tt(7}-rM) 

h\ 

(41) 

ui+'/i(T  +  ~  ui-y,(T  +  Nl)/^ 
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wi+  h(T  +  N1  )/+  Vl  -  (T  +  Nl )j_Vx 

K 

(T  +  N 1 )/+  ft ,  etc.,  are  computed  using  upwind  differencing  as 
with  Ti+Vl  in  eq  39.  iV2(+Vi,  etc.,  are  defined  in  a  manner  simi¬ 
lar  to  the  velocities  defined  in  eq  38,  i.e.,  the  average  of  the  two 
nodes  adjacent  to  the  nodal  boundary  is  taken.  Velocities  at 
the  node  boundaries  are  computed  as  in  eq  38. 

Applying  the  finite  difference  formulae  derived  from  Taylor 
series  is  much  more  difficult  at  the  boundaries,  and  thus  the 
control  volume  approach  proves  superior  there.  Next  we  shall 
illustrate  the  approach  for  an  adiabatic,  impermeable  side  wall. 
Complete  derivations  for  all  of  the  cases  of  interest  are  given  in 
Appendix  A. 

The  control  volume  is  defined  in  Figure  8 

The  boundary  condition  is  that  no  h  t  •  s  through  the 
wall.  Thus  the  compc  mts  of  the  ener  ition  are 
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Figure  8.  Control  volume. 


The  control  volume  used  for  interior  nodes  is  illustrated  in  Figure  7.  Earlier  we  derived  an 
energy  equation  for  single-component  fluid  flowing  through  a  porous  medium.  The  physical 
statement  of  eq  30  is  that 

net  (conduction  +  convection)  =  accumulation. 

It  was  assumed  that  the  system  was  at  steady  state,  so  that  accumulation  was  zero.  Here  we  in¬ 
troduce  a  transient  term  so  that  we  may  numerically  advance  the  system  to  its  final  steady  state 
For  the  control  volume  about  the  node  ij  the  accumulation  over  one  time  step  At  is  written  as 

acc  =  (T*f 1  -  Tfj )hxhz 


where  the  index  C  indicates  the  time  step  level,  and  T,  x  and  z  are  the  nondimensional  variables 
defined  earlier. 

The  flow  of  heat  into  the  control  volume  by  conduction  in  the  time  At  is  given  by 


Physically,  the  conduction  terms  are  positive  when  the  temperature  is  greater  outside  the  control 
volume  than  inside,  i.e.,  heat  flows  in.  The  heat  flow  across  each  boundary  is  simply  the  gradient 
across  that  boundary,  which  here  is  approximated  with  a  central  difference  applied  at  the  con¬ 
trol  volume  boundary,  multiplied  by  the  area  for  heat  transfer  and  the  time  interval. 

For  the  convective  terms,  heat  flows  in  when  the  fluid  flows  in.  The  heat  flow  attributable  to 
convection  is 


convect  =  ((u D,_w  -  <«D(4V)  |h.Ar  -  (<*vD;_Vl  -(*'7’)/+,/,  ]fcx At 
where  u  and  w  are  positive  in  the  positive  x  and  z  directions  respectively.  We  define 


V/-i  +  *../ 


(38) 


The  upwind  differencing  concept  then  defines 


for  ut  s  >  0 

T  . 

1  i./-i 

for  >  0 

= 

- 

J,, 

for  <  0 

J,, 

for  <  0 

(39) 

T 

for  <  0 
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for  w/>%  >  0 

T,+  V, 
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Tj+'/l  = 
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o 

7u*  i 

for  w/+ 1/4  <  0  . 

In  calculating  the  heat  fluxes  across  each  boundary,  we  have  assumed  that  the  temperature 
and  temperature  gradient  normal  to  the  control  volume  boundary  are  constant  along  the  bound- 
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because  the  value  for  i jj  is  known.  However,  when  the  top  is  permeable,  the  boundary  condition 
is  d\p/dz  =  0,  and  we  must  solve  for  \p  at  the  upper  boundary.  We  cannot  simply  apply  eq  34  be¬ 
cause  the  node  at  index ;+l  does  not  exist  for;  =  nz.  The  procedure  used  is  to  employ  a  shadow 
node  technique,  in  effect  pretending  that  the  node  is  there.  We  write  a  central  difference  form 
of  the  boundary  condition  as 


^i,nz  +  1  ~  ^i.nz-l 

2  hz 


=  0 


(35a) 


and  a  backward  difference  form  of  the  temperature  derivative  as 

(df\  3Ti,nz  -4Ti,nz- 1  +  Ti,nz- 2 

2h, 


(35b) 


Then,  by  writing  the  PDE  at  j=nz,  we  can  eliminate  \l/fl2+ ,  by  solving  for  it  from  eq  35a.  The  re¬ 
sult  is  an  equation  for  i pnz : 

~  2\j/;  'Pj-i  )nz  ^  ^  nz-1  ~  ^i,nz  _ 


Ra 


r  ,(^i+  1  ,nz  ~  ^i-l,nz\  .  '47'-_|  +  Tt  _^\ 

lc0*V — nr, — )  ■  ”*( - nr, - j 


(36) 


Velocities  are  calculated,  again  using  three-point  central  differences.  For  interior  nodes 


hi*  t  -  1 


1.7 


2  h. 


(37) 


^>1,7 1,7 


WU  =  -  2  h 


At  the  boundaries,  the  velocity  component  normal  to  the  wall  is  zero,  but  the  component  parallel 
to  the  wall  is  not.  The  appropriate  three-point  backward  or  forward  difference  is  then  employed, 
similar  to  eq  35b. 

To  derive  the  discrete  form  of  the  energy  equation,  integration  is  performed  over  a  control  vol¬ 
ume  about  each  node.  The  discrete  form  may  also  be  found 
from  Taylor  series  formulas,  but  this  method  is  physically 
less  appealing  and  difficult  to  apply  at  the  boundaries.  Up¬ 
wind  differencing  is  applied  to  evaluate  the  convection  terms. 

Essentially,  this  means  that  only  nodes  upstream  of  the  con¬ 
trol  volume  boundary  are  allowed  to  influence  the  heat  trans¬ 
fer  across  the  boundary.  Upwind  differencing  is  used  to 
damp  the  instability  associated  with  the  non-linear,  first  or¬ 
der  convective  terms.  The  method  enhances  the  stability  of 
the  numerical  method,  but  at  a  cost  in  accuracy,  because  the 
upwind  difference  expressions  have  a  truncation  error  of 
CH h).  This  concept  will  be  illustrated  shortly. 


-  4 - 


-t.)  ,  ;.i 

I 

1  >,i-i 


i»t.i 


Figure  7.  Control  volume  used 
for  interior  nodes. 
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From  continuity  for  the  vapor  phase,  rv  =  V'CPv^)  +  V‘/v  at  steady  state,  and  by  definition, 

L  =  hv  -  hs.  The  gradients  of  pv  as  a  function  of  temperature  are  obtained  by  differentiating  the 
equation  of  state.  Using  these  relations,  and  Fick’s  law  for  jv,  we  find  the  final  dimensional  form 
of  the  energy  equation 

0 0Cp\  V  •  vT  +  L  V  -pv?  -  2  T  +  LBD  V  •  (pvVT).  (30) 

Finally,  we  put  eq  26  and  eq  30  in  nondimensional  form.  T,  V  and  v  are  the  dimensional  quan¬ 
tities  that  are  now  used  in  nondimensional  form,  which  is  done  identically  to  that  in  the  Equation 
of  Motion  section.  Then 

V’vT  =  V2  T  (31) 


and 


^•(r  +  JVOi^V  \N2VT) 
where 


(32) 


N\ 


(f>Cp)aAT 


and  N2  =  1  +  • 


LDBpv 


It  is  important  to  note  that  both  N\  and  N2  are  functions  of  temperature. 


(33) 


Finite  difference  methods 

Finite  differences  were  chosen  as  a  means  of  obtaining  approximate  solutions  to  the  partial 
differential  equations  (PDE).  The  method  breaks  the  region  of  interest  down  into  discrete  points, 
replacing  the  PDE  with  a  system  of  algebraic  equations  in  terms  of  values  of  the  dependent  vari¬ 
ables  at  the  discrete  points.  As  the  number  of  points  increases,  both  the  accuracy  of  the  solution 
and  the  computing  time  increase.  Finite  differences  are  especially  well  suited  to  the  present  prob¬ 
lem  because  of  the  simple  geometries  involved.  A  typical  grid  for  the  current  problem  is  illustrat¬ 
ed  in  Figure  6. 

Formulas  for  discretizing  the  PDEs  can  be  derived  either  from  Taylor  series  expansions  (Pinder 
and  Gray  1977)  or  by  integrating  over  control  volumes  around  nodes.  This  latter  method  is  espe¬ 
cially  useful  for  conservation  equations  such  as  the  energy  equation  in  the  current  problem.  We 
later  demonstrate  this  approach,  which  is  used  extensively  in  deriving  formulae  for  the  energy 
equation. 

First,  we  will  present  the  results  for  the  equation  of  motion.  Finite  difference  expressions, 
which  have  a  truncation  error  of  0(h2 ),  are  used  for  all  terms;  h  is  the  spacing  between  grid  points. 
The  discrete  form  of  eq  14  is 


1 ,/  -  2*u +  *<-!,/ ,  *<,/+ 1  -  2^U  +  *u- 1 
hi  hi 


(34) 


This  form  is  valid  at  all  interior  nodes. 

When  the  boundary  condition  is  the  Dirichlet  condition 
(< p  =  0),  we  do  not  need  to  apply  the  PDE  at  the  boundary 


j  *nz 


i  =  l 


i  =  nx 


Figure  6.  Typical  finite  dif¬ 
ference  grid. 


We  can  evaluate  the  enthalpy  by  the  following  thermodynamic  relations: 


dh  =(i)  dP  •  m  =cp  and  ®  -en 


V>P It  p 


where  <3  is  the  isobaric  coefficient  of  thermal  expansion.  For  air  0  *=  1/7',  and  thus 


dh  =  CpdT .  (24) 

Applying  Fourier’s  law  of  heat  conduction: 

q  =  -km$T  (25) 

and  the  energy  equation  is,  in  final  form, 

(pCp)ai?’vT=  kmV2T .  (26) 

This  is  the  standard  form  of  the  energy  equation  used  in  modeling  of  porous  media. 

We  now  introduce  vapor  effects  so  pv  varies  with  temperature  and  r{  is  not  zero.  Applying 
these  to  eq  22,  along  with  rg  =  -rv,  results  in 

/VVVM  Pava4ha  +rv(hv-hs)  =  km?iT.  (27) 

At  steady  state  rv  =  A(pv)v.  Component  velocities  are  given  by 
vv  =  v+!jpv 
\  ?  +LlPa  - 

The/’s  are  the  component  velocities  relative  to  the  mass  average  velocity.  Here  they  are  consid¬ 
ered  to  be  caused  by  concentration  diffusion  only.  Since  ja  =  -  jv,  and  from  Fick’s  law  of  diffu¬ 
sion,  jv  =  -  DVPv,  we  can  estimate  the  relative  contributions  of  diffusion  to  the  component  veloc¬ 
ities.  For  an  equation  of  state  for  the  vapor,  we  use  the  exponential  relation 

pv=p0exp\B(T' -Ta)\  .  (28) 

Here  we  have  assumed  that  the  vapor  is  everywhere  saturated.  Although  this  is  not  strictly  true, 
typical  supersaturations  are  usually  less  than  0.1%  and  thus  have  little  effect  on  the  overall  diffu¬ 
sion  process. 

At  0°C,  with  a  0.3°C/cm  temperature  gradient,  the  diffusional  flux  is  about  2.7x  10-8  g/cmJ  s. 
The  diffusional  velocities  are  then  5x  10"3  and  2x  10“*cm/s  for  vapor  and  air  respectively.  Since 
convective  velocities  are  of  the  order  of  0.01  cm/s,  it  is  certainly  safe  to  neglect  at  least  the  diffu¬ 
sional  component  of  air.  Also,  the  density  of  vapor  is  three  orders  of  magnitude  less  than  that  of 
air,  and  since  the  heat  capacities  are  similar,  the  second  term  of  eq  27  is  much  greater  than  the 
first.  The  energy  equation  is  now  of  the  form 

(pCp)aV-vT+rv(hv-h3)=kmV2T.  (29) 


M 

r-v.-q 
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term.  The  gravitation  term  is  estimated  to  be  four  orders  of  magnitude  less  than  the  comparable 
convection  term,  and  thus  is  also  neglected. 

Since 

=  A,  -  (Plp)i 
we  can  write  eq  16  as 

S^(pA-p),+  v-/™]/='V‘9-V‘P>'  (,7 

where  the  summation  is  taken  over  all  components,  which  may  have  different  velocities  because 
they  diffuse  at  different  rates.  Since 

dp 

Pi  =  P,/2p,-,  2  V  •  (pv)t  =  V  •  pv  and  2  —  = 
the  two  pressure  work  terms  are  identical  and  cancel  each  other  out.  Rearranging,  we  get 


?  hi(p  +v-  p?), +  ?(<>  ff V  =  -V*  q  ff  • 

The  continuity  equation  for  each  component  is 


=  -V-(P^  \+r‘ 


where  r(  is  the  rate  of  production  of  component  i.  The  equation  of  continuity  for  an  incompres¬ 
sible  fluid  is 

V*  v  =  0.  (20 


Combining  eq  1 8  and  19  results  in 


(  „  ,\  dp 


Here  we  deal  only  in  the  steady-state  results,  assuming  that  the  phenomena  in  snow  are  at  least 
quasi-steady,  that  is,  changing  slowly. 


This  is  the  general  form  of  the  energy  equation  we  will  apply. 

First  let  us  consider  only  air  flowing  through  snow,  with  no  vapor  contribution  (or  vapor  con 
tributions  independent  of  temperature).  Then  r,-  =  0  and  v  =  v  The  resulting  energy  equation 


If  we  take  the  curl  (Vx)  of  both  sides,  the  pressure  term  is  identically  zero  and  the  result  is 


3w 

3* 


3  u  KH 

5 


where  w  and  u  are  the  dimensionless  velocity  components  in  the  z  and  x  directions  respectively. 

To  a  good  approximation,  the  density  of  air  is  a  linear  function  of  temperature  over  the  tem¬ 
perature  range  of  interest,  thus  we  can  write  an  equation  of  state  of  the  form 


Pf=PoV-(Kr  -  t0)]. 


(10) 


By  introducing  a  dimensionless  temperature  defined  as 


T'  -T 
T  -  o 

AT 


(ID 


and  substituting  this  into  the  density  term,  the  equation  of  motion  becomes 


3w  3  m 

____= - - - (VxgT) . 

3x  3z  hk 

Finally,  a  stream  function  of  the  form 


(12) 


dil/  dip 

—r=w  —X.  -  -u 

3x  3z 


(13) 


is  introduced.  The  final  form  of  the  equation  of  motion  is  then 


dT  bT 
V 2  Ip  =  Ra  [cos0  ^  -  Sin0  ] 


(14) 


where 


Pog0ATKH 
Ra  = - — — 

fiK 

As  discussed  earlier,  the  parameter  AT  is  defined  by  the  boundary  conditions. 

Energy  equation 

For  a  multicomponent  system,  the  energy  equation  can  be  written  as 


S  §f(pV)t  +  2  V  •  (P ?U)t  =  V  •  q -V  •  pt  +  Pp{ v  ■  g\-  V •  ( t- •  v) 


(15) 


(16) 


where  the  summation  is  over  all  components.  Here  primes  are  not  used,  although  the  variables 
are  later  made  nondimensional.  This  is  essentially  a  multicomponent  form  of  an  equation  pre¬ 
sented  by  Bird  et  al.  (1960).  All  of  the  work  terms  have  been  written  for  a  single-component 
fluid,  because  for  the  snow-air  system  the  fluid  is  dominantly  air.  Kinetic  energy  terms  have  been 
neglected.  The  last  term  on  the  right  side  represents  viscous  work  ascribable  to  energy  dissipation 
in  the  fluid  but,  since  dissipation  is  dominated  by  fluid  shear  at  the  solid  matrix,  we  neglect  this 


accounting  for  latent  heat  effects  from  vapor  transport.  Both  models  assume  that  the  medium 
is  isotropic  and  homogeneous.  Boundary  conditions,  aspect  ratio  of  the  medium,  Rayleigh  num¬ 
ber  and  slope  may  all  be  varied  for  both  models.  In  addition,  for  the  second  model  the  effective 
vapor  diffusivity  (Deff)  and  thermal  conductivity  (km )  may  also  be  varied. 

In  this  section,  the  governing  equations  are  developed  and  the  numerical  methods  are  outlined; 
in  the  next  section  the  results  are  presented. 

Equation  of  motion 

It  is  assumed  here  that  air  flow  through  snow  can  be  described  by  the  empirical  relation  devel¬ 
oped  by  Darcy.  This  is  valid  for  Reynolds  numbers  (Re)  less  than  a  critical  number  for  the  onset 
of  turbulent  flow,  which  falls  between  1 .0  and  1 0.0  (Dullien  1979).  The  Reynolds  number  is  de¬ 
fined  as 


Re  = 


Pfd 

77 


where  d  is  a  particle  diameter.  Our  calculations  (see  Applications  to  Snow  Metamorphism  section) 
reveal  that  a  typical  velocity  is  0.04  cm/s.  Thus  for  a  1-mm  particle,  the  Reynolds  number  is 
about  0.03,  and  Darcy’s  Law  should  be  valid. 

Darcy’s  Law  is  traditionally  (Scheidegger  1974)  written  as 

V'p- |  v' -  pfg  =  0  (8) 

where  g  is  the  gravitational  vector  and  is  in  the  same  direction  as  the  positive  z  axis.  The  velocity 
is  averaged  over  a  cross-sectional  area  that  is  an  order  larger  than  the  pore  size  scale,  thus  the  ve¬ 
locity  is  equivalent  to  a  volumetric  flux.  When  the  fluid  is  i  .ade  up  of  several  components,  a 
mass  average  velocity  is  used,  defined  as 


?(p*0, 

r 


where  the  summation  is  over  all  fluid  components  and  each  v(-  is  the  volumetric  flux  of  the  compo¬ 
nent.  The  primes  indicate  quantities  later  made  dimensionless. 

We  limit  ourselves  to  two-dimensional  modeling  and  define  the  coordinate  system  in  two  di¬ 
mensions  as  shown  in  Figure  5.  The  gravity  vector  has  negative  z  and  x  components  and  thus 
the  final  term  of  eq  8  has  the  opposite  sign.  Introducing 
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(9) 


Figure  5.  Two-di¬ 
mensional  coordi¬ 
nate  system. 


Studies  of  convection  through  snow 

To  date,  only  Akitaya  (1974)  and  Palm  and  Tveitereid  (1979)  have  studied  thermal  convection 
in  snow.  Akitaya ’s  work  was  primarily  experimental,  while  that  of  Palm  and  Tveitereid  was  en¬ 
tirely  theoretical. 

Akitaya  attempted  to  generate  convection  in  a  layer  of  snow  15  cm  deep,  but  was  largely  un¬ 
successful.  He  was  able  to  generate  convection  only  in  artificial  snow  with  an  average  grain  size 
of  1 5  mm.  That  he  was  able  to  generate  convection  in  only  snow  of  extreme  permeability  is  not 
surprising,  since  convection  is  less  likely  in  shallow  snow  layers.  His  results  with  artificial  snow 
are  surprising  because  they  show  a  step  change  in  the  effective  thermal  conductivity  with  increas¬ 
ing  temperature  gradient,  rather  than  the  continuous  increase  we  would  expect. 

Akitaya ’s  comparison  of  experiment  with  theory  led  him  to  conclude  that  in  snow,  convection 
occurs  at  a  much  higher  than  predicted  critical  gradient.  His  conclusion  is  not  valid,  however,  be¬ 
cause  he  used  pCp  for  snow  rather  than  for  air  when  he  calculated  the  conditions  for  the  onset  of 
convection.  When  correctly  calculated,  it  is  obvious  why  there  was  no  convection  in  his  natural 
snow  experiments. 

Palm  and  Tveitereid  (1979)  calculated  the  critical  conditions  for  convection  in  snow  in  a  simi¬ 
lar  manner  to  that  used  for  other  porous  media.  In  general,  their  conclusion  was  that  convection 
could  occur  only  in  old,  coarse-grained  snow  subjected  to  severe  gradients.  While  it  is  true  that 
these  are  conditions  under  which  convection  is  most  likely,  they  probably  underestimate  the  like¬ 
lihood  of  convection  under  less  severe  situations.  In  particular,  they  used  a  value  of  air  permea¬ 
bility  for  “coarse-grained”  snow  that  is  similar  to  that  observed  by  Shimizu  (1970)  for  “fine¬ 
grained”  snow.  Shimizu  reported  values  about  twice  as  large  for  what  he  calls  coarse-grained  snow, 
and  Bader  (1939)  observed  permeabilities  for  depth  hoar  as  much  as  five  times  the  largest  consid¬ 
ered  by  Palm  and  Tveitereid. 

The  theory  developed  in  the  sections  above  is  not  sufficient  to  fully  describe  convection  in 
snow.  Phase  change  from  the  net  condensation  or  sublimation  of  water  vapor  is  an  important 
part  of  the  heat  transfer  process.  This  was  recognized  by  Palm  and  Tveitereid,  but  was  incorporat¬ 
ed  into  their  analysis  only  insofar  as  the  value  of  k  was  affected  in  the  Rayleigh  number.  Since 
their  value  of  km  was  arbitrarily  chosen,  their  consideration  of  phase  change  was  physically  mean¬ 
ingless.  It  is  necessary  to  treat  phase  change  as  a  part  of  the  governing  equations,  since  its  influ¬ 
ence  depends  on  local  temperatures,  temperature  gradients  and  velocities. 

The  boundary  conditions  in  snow  frequently  include  a  permeable  top  and  a  constant  heat  flux 
bottom.  While  the  influence  of  these  boundary  conditions  on  the  onset  of  convection  is  well  un¬ 
derstood  (see  Table  1 ),  their  influence  on  the  intensity  of  convection  at  Rayleigh  numbers  greater 
than  critical  has  not  been  quantified.  This  is  one  of  the  aims  of  later  sections  of  this  report. 

Those  recently  working  on  convection  in  snow  have  ignored  the  convection  phenomenon  asso¬ 
ciated  with  sloped  layers.  Neher  (1939)  recognized  that  flow  of  air  occurs  up  a  snow  slope  when 
a  temperature  gradient  is  applied,  but  none  of  his  successors  have  made  mention  of  this  fact.  He 
attributes  the  existence  of  thicker  depth  hoar  layers  in  the  upslope  parts  of  inclined  snow  covers 
to  this  phenomenon.  This  observation  would  seem  of  vital  importance  in  avalanche  predictions 
because  these  upslope  regions  are  frequently  starting  zones  for  avalanches. 

In  the  following  section,  we  address  the  above-stated  gaps  in  our  current  knowledge.  Later  our 
findings  are  applied  to  two  real  problems  of  convection  in  a  snow  layer:  First,  does  it  occur?  And 
second,  is  it  important  in  snow  metainorphism? 


MODELING 

A  numerical  model  is  developed  here  to  simulate  convection  in  a  porous  medium.  Two  forms 
of  the  model  are  developed,  one  form  without  phase  change,  which  is  similar  to  the  porous  media 
model  of  Ribando  (1977).  The  second  form  more  closely  models  convection  through  snow  by 


Figure  4.  Unicellular  convection  in  a  sloped  porous  layer  (after  Combar- 
nous  and  Bories  1975). 

fer  is  affected.  This  convection  is  illustrated  in  Figure  4;  the  velocity  profile  can  be  described  by 
(Bories  and  Combarnous  1973) 

u(z)  =  PoPg  A  TK  sin  -j^j  .  (7) 

Thus  the  intensity  of  convection  increases  with  increasing  slope. 

As  the  Rayleigh  number  increases  to  the  critical  value,  the  structure  of  convection  becomes 
more  complex.  For  slightly  inclined  slopes  (<t>  <  10-1 5°),  hexagonal  cells  have  been  observed  by 
Bories  and  Combarnous  (1973).  Kaneko  et  al.  (1972)  have  observed  transverse  rolls  (axis  parallel 
to  the  contour)  at  similar  low  angles.  The  reason  they  observed  rolls  and  not  hexagonal  cells  was 
that  the  dimensions  of  their  box  were  much  smaller  along  the  roll  axis,  and  thus  two  dimensional 
convection  was  forced. 

In  more  steeply  inclined  layers  at  high  Rayleigh  numbers,  Bories  and  Combarnous  found  long¬ 
itudinal  rolls  (axis  parallel  to  fall  line).  Kaneko  et  al.  observed  a  single  unicellular  pattern  in  steep 
layers.  The  two  sets  of  observations  are  consistent  when  respective  container  geometries  are  con¬ 
sidered.  Krantz  et  al.  (1983)  observed  a  similar  transition  in  field  studies  of  patterned  ground. 
The  transition  they  observed  takes  place  between  3°  and  7°. 

The  behavior  of  these  cells  (or  rolls)  is  similar  to  that  in  horizontal  systems,  except  that  the 
governing  parameter,  Ra,  is  replaced  by  Ra  costp.  The  experimental  data  of  Bories  and  Combar¬ 
nous  (1973)  show  that  if  the  Nusselt  number  is  plotted  as  a  function  of  Ra  c os<p,  the  result  is  iden 
tical  to  Figure  3,  especially  at  high  Rayleigh  numbers.  Weber  (1975)  extends  the  analytical  re¬ 
sults  of  Palm  et  al.  (1972)  by  simply  replacing#  by  g  cos^i  in  the  definition  of  Rayleigh  number. 
The  analytical  results  agree  well  with  the  experimental  data  of  Dories  and  Combarnous. 

No  work  has  been  done  on  the  influence  of  different  boundary  conditions  or  the  effects  of 
lateral  confinement  on  convection  in  sloped  layers.  Also,  the  existence  of  the  longitudinal  rolls 
described  by  Bories  and  Combarnous  (1973),  superimposed  on  the  basic  uphill-downhill  unicell¬ 
ular  flow,  must  be  better  understood.  No  return  path  for  the  fluid  is  described  and  thus  the  spi¬ 
ral  coil  may  not  exist. 


Table  2.  Summary  of  the  results  of  Ribando  (1977),  indicating  the 
influence  of  constant  flux  bottom  and  permeable  top  boundary 
conditions. 


Boundary  conditions 

1 Top 

Bottom  AR 

Ra  Nu 

^max 

Constant  temperature, 

Constant  temperature,  0.6 

20C  4.07 

6.87 

impermeable 
Constant  temperature. 

impermeable 

Constant  temperature,  0.6 

2  6.38 

9.06 

permeable 

Constant  temperature. 

impermeable 

Constant  flux,  0.6 

200  =  2.6 

3.1 

impermeable 
Constant  temperature. 

impermeable 

Constant  flux,  0.6 

200  «  3.1 

3.21 

permeable 

Constant  temoerature. 

impermeable 

Constant  flux,  0.6 

3000  NA* 

7.98 

impermeable 
Constant  temperature. 

impermeable 

Constant  flux,  0.6 

3000  NA 

13.85 

permeable 

Constant  temperature. 

impermeable 

Constant  flux,  1.0 

H 

o 

TT 

1.44 

permeable 

Constant  temperature. 

impermeable 

Constant  flux,  1.0 

100  ^  1.9 

2.86 

rrmeable _  impermeable 


♦Not  available 

ability  of  the  medium  to  transfer  heat  has  been  doubled.  When  the  boundary  conditions  are 
fixed  temperatures,  the  result  is  a  doubling  of  the  average  heat  flux.  When  one  of  the  boundaries 
is  a  fixed  flux  boundary,  then  the  average  temperature  difference  necessary  to  maintain  a  given 
heat  flux  is  halved.  Figure  3  shows  some  of  the  results  to  date,  both  experimental  and  theoreti¬ 
cal,  for  heat  transfer  between  isothermal,  impermeable  boundaries. 

Experimental  results  for  other  boundary  conditions  were  not  obtained.  Ribando  (1977)  in¬ 
vestigated  numerically  the  effects  of  a  permeable  top  and  of  a  constant  flux  bottom  upon  heat 
nsfer  in  porous  media  (Table  2).  In  general,  his  results  indicate  that  for  a  given  Rayleigh  num¬ 
ber  a  permeable  top  increases  the  intensity  of  convection,  while  a  constant  flux  bottom  decreases 
it,  compared  to  isothermal,  impermeable  boundaries.  It  is  apparent  that  the  correct  choice  of 
boundary  conditions  is  very  important,  yet  doing  so  is  often  difficult,  especially  since  constant 
flux  and  isothermal  boundaries  are  identical  at  subcritical  Rayleigh  numbers. 

Layering  and  slope  effects 

The  effects  of  layering  and  slope  are  of  great  importance  to  this  study,  since  snow  covers  are 
frequently  heterogeneous  and  are  often  found  on  hillsides.  Most  of  the  analysis  of  the  effects  of 
layering  has  been  done  using  numerical  solutions  of  the  governing  differential  equations.  Both 
the  onset  of  convection  (McKibben  and  O’Sullivan  1980,  Richard  and  Gounot  1981)  and  heat 
transfer  in  layered  porous  media  (Rana  et  aJ.  1979,  McKibben  and  O’Sullivan  1981 ,  Richard  and 
Gounot  1981)  have  been  treated.  Some  of  the  results  have  been  compared  favorably  with  exper¬ 
iments  (Richard  and  Gounot  1981)  and  with  data  from  geothermal  reservoirs  (Rana  et  al.  1979). 

Because  of  the  infinite  number  of  possibilities  for  combinations  of  layers,  generalized  correla¬ 
tions  are  not  available  to  predict  the  effects  of  layers  on  the  onset  of  convection  and  heat  trans¬ 
fer  by  convection.  McKibben  and  O’Sullivan  (1980)  conclude  that  unless  permeability  differences 
between  layers  are  50%  or  more,  the  onset  of  convection  is  not  affected.  Also,  when  an  upper 
layer  of  low  permeability  exists,  a  closed  or  open  upper  boundary  does  not  affect  the  results. 

The  outstanding  feature  that  distinguishes  the  situation  in  a  sloped  layer  from  that  in  a  hori¬ 
zontal  layer  is  that  on  a  slope  there  will  be  some  motion  of  air  for  any  finite  A T.  At  subcritical 
Rayleigh  numbers,  the  flow  is  a  simple  unicellular  pattern,  rising  along  the  lower  warm  boundary 
and  returning  downward  below  the  cold  upper  boundary.  In  a  slope  of  infinite  extent,  flow  per¬ 
pendicular  to  the  slope  is  negligible,  and  thus  neither  the  temperature  profile  nor  the  heat  trans- 


The  critical  Rayleigh  number  is  minimized  for  a  given  set  of  boundary  conditions  at  some  op¬ 
timum  cell  size.  The  critical  wave  number  acr  is  also  given  in  Table  1 .  This  wave  number  is  de¬ 
fined  as 

alr  =  %2+m2 

where  £  and  m  are  the  wave  numbers  in  the  x  and  y  directions,  defined  respectively  as  £  =  2 nH/Wx 
and  m  =  2nH/Wy ;  Wx  and  Wy,  are  the  cell  dimensions  in  the  x  and y  directions  respectively. 
Knowing  acr,  we  can  predict  the  cell  size  at  the  onset  of  convection  for  either  two-dimensional 
rolls  or  three-dimensional  hexagons. 

The  above  results  apply  when  the  lateral  extent  of  the  layer  is  large.  When  confinement  is  im¬ 
portant,  the  dimensions  of  the  container  dictate  the  cell  form  and  the  critical  Rayleigh  number. 
Beck  (1972)  and  Tewari  and  Torrence  (1981)  solved  for  the  cell  form  and  critical  Rayleigh  num¬ 
ber  as  a  function  of  container  size  for  the  closed  and  open  top  cases  respectively.  They  found 
that  unless  both  the  container  dimensions  are  much  less  than  the  critical  wavelengths,  the  critical 
Rayleigh  number  is  not  increased  significantly.  We  aiso  conclude  that  when  one  of  the  lateral 
dimensions  is  less  than  the  critical  wavelength,  the  preferred  cell  form  at  the  onset  is  two-dimen¬ 
sional  rolls. 

Heat  transfer  attributable  to  thermal  convection 

When  the  Rayleigh  number  increases  beyond  the  critical  value,  convection  occurs,  causing  a 
significant  increase  in  the  rate  of  heat  transfer.  Numerical,  analytical  and  experimental  methods 
have  all  been  used  to  investigate  this  phenomenon.  Published  results  are  commonly  presented  in 
terms  of  the  Nusselt  number,  defined  as 

Nu  =  -&■  (6) 

where  keff  is  an  effective  thermal  conductivity  that  accounts  for  both  conductive  and  convective 
heat  transfer.  If  Nu  =  1 .0,  convection  has  no  effect  on  heat  transfer,  whereas  if  Nu  =  2.0,  the 


Figure  3.  Experimental  and  theoretical  results  for  heat  transfer  bet  wen 
isothermal,  impermeable  boundaries. 
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which  is  valid  for  Ra/Racr  >  1 ,  is  shown 
for  comparison.  We  conclude  that  the  re¬ 
lation  of  Elder  describes  the  experimental 
results  reasonably  well.  Figure  10  shows 
numerical  results  from  our  model  for  the 
same  boundary  conditions  as  Figure  9, 
and  makes  a  favorable  comparison  with 
eq  43.  The  runs  leading  to  these  results 
were  done  with  AR  =  1 ,  which  is  equiva¬ 
lent  to  the  cell  size  at  the  onset  of  con¬ 
vection.  Runs  are  made  only  for  Ray¬ 
leigh  numbers  just  above  the  critical  val¬ 
ue,  where  the  cell  size  is  not  too  differ-  o  i  2  3 

ent  from  that  at  the  onset.  By  compar-  RO/Rocr 

ing  the  results  of  Elder  to  the  experi¬ 
mental  data,  we  conclude  that  our  nu-  Figure  10.  Numerical  results  and  the  relation  de- 
merical  model,  which  agrees  with  eq  43,  duced  by  Elder  ( 1 967)  for  convection  between  iso- 
is  actually  a  better  tool  for  predicting  thermal,  impermeable  boundaries. 

heat  transfer  rates  than  either  of  the 
analytical  results. 

When  the  bottom  boundary  is  a  constant  flux  thermal  condition,  we  have  only  the  numerical 
results  of  Ribando  (1977)  for  comparison.  Table  3  shows  numerical  results  from  our  model  with 
those  of  Ribando  for  identical  conditions.  Nusselt  numbers  for  Ribando  are  estimated  from  his 
plots  of  isotherms,  using  the  temperature  difference  at  the  midpoint  of  the  top  and  bottom  bound 
aries.  Nusselt  numbers  for  the  present  model  are  calculated  by  integration  along  the  boundaries 
to  calculate  an  average  AT.  Agreement  for  all  cases  is  very  good. 

Finally,  we  test  the  model  by  trying  to  predict  the  onset  of  convection.  This  is  done  for  the 
closed  top  isothermal  and  open  top  isothermal  cases.  Aspect  ratios  used  are  those  of  a  single  cell 
size  at  the  onset  of  convection,  calculated  from  the  critical  wavelength  values  given  in  Table  1. 
While  it  is  not  really  practical  to  identify  an  exact  onset  point,  it  is  possible  to  find  a  narrow  range 
in  which  the  onset  point  lies.  The  upper  and  lower  limits  indicate  Rayleigh  numbers  at  which 
convection  clearly  was  and  was  not  occurring  respectively.  We  find  that  38.0  <  Racr  <  41 .0  and 
26.0  <  Racr  <  28.0  for  the  closed  and  open  top  cases  respectively.  The  corresponding  analytical 
results  of  Ra  =  39.5  and  Ra  =  27.1  fit  within  the  observed  range  for  each  case. 

The  effects  of  phase  change  have  yet  to  be  dealt  with  in  the  literature,  and  thus  no  means  of 
testing  the  model  for  significant  phase  change  exist.  Experimental  and  theoretical  results  exist 
for  convection  in  a  sloped  layer,  and  the  model  predicts  the  general  trends  indicated  by  the  pub- 

Table  3.  Comparison  of  published  values  with  this  study. 

All  values  are  based  on  numerical  calculations.  Nusselt 
numbers  for  the  results  of  Ribando  are  estimated  from 
plots  of  isotherms.  (After  Ribando  [1977J  and  Klever 
[1983].) 
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Open  top,  constant 

1.0 
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1.44 

Ribando 

flux  bottom 

1.0 

40 

1.41 1 

1.44 

Present  work 

Closed  top,  constant 

1.0 

100 

*1.9 

2.86 

Ribando 

flux  bottom 

1.0 

100 

1.951 

2,82 

Present  work 

Closed  top. 

0.6 

200 

4.07 

6.87 

Ribando 

isothermal  bottom 

0.6 

200 

NA 

6.74 

Klever 

0.6 

200 

4.02 

6.69 

Present  work 
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lished  works.  However,  the  manner 
of  presentation  of  results  in  the  lit¬ 
erature  obscures  some  of  the  physi¬ 
cal  phenomena  that  occur  in  sloped 
layers.  We  have  shown  here  that  the 
model  is  a  reliable  tool  for  analyz¬ 
ing  convection  in  horizontal  porous 
media  and  it  should  allow  us  to  study 
cases  for  which  little  data  are  current¬ 
ly  available. 


a.  Closed  top,  isothermal  bottom. 


MODELING  RESULTS 

The  effects  of  boundary  condi¬ 
tions  on  heat  transfer  and  air  veloc¬ 
ity  are  studied  using  the  model. 
Subsequently,  the  effects  on  heat 
transfer  attributable  to  latent  heat 
release  by  phase  change  are  analyzed 
and  found  to  be  a  function  of  the 
diffusivity.  Finally,  we  use  the  mod¬ 
el  to  investigate  several  issues  of  im¬ 
portance  in  sloped  layers,  including 
the  effects  of  lateral  confinement. 


X=l.35 


b.  Open  top,  isothermal  bottom. 


Effects  of  constant  flux  and 
permeable  boundaries  on 
convection  in  horizontal  layers 

The  work  of  Ribando  (1977)  is 
an  important  starting  point  in  the 
discussion  of  the  effects  of  boundary 
conditions  upon  convection  in  hori¬ 
zontal  layers.  He  developed  the 
mathematical  forms  of  the  boundary 
conditions  of  interest,  and  applied 
them  to  several  simple  cases.  The 
aim  here  is  to  more  thoroughly  in¬ 
vestigate  the  effects  of  these  boun¬ 
dary  conditions,  so  that  the  more 
general  conclusions  about  the  effects 
on  heat  transfer  can  be  made. 

Figure  1 1  shows  the  steady-state 
results  of  runs  for  Ra/Racr  =  1 .5  for 
four  combinations  of  boundary  con¬ 
ditions.  The  aspect  ratio  in  each  case 
corresponds  to  the  size  of  a  single 
cell  at  the  onset  of  convection,  for 
the  given  boundary  conditions.  When 
the  Rayleigh  number  is  greater  than 
the  critical  value,  the  size  of  a  convec- 


X  =  I35 


c.  Closed  top,  flux  bottom. 


— iz=o 

X  =  l.8 


d.  Open  top,  flux  bottom. 

Figure  11.  Isotherms  (solid  lines )  and  streamlines 
(dashed  lines)  for  Ra/Racr  =  1-5  for  boundary 
conditions  as  labeled. 
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Figure  12.  Heat  transfer  re¬ 
sults  of  the  numerical  mod¬ 
eling  in  this  report,  and  some 
results  of  Elder  (1967)  and 
Ribando  (1977). 


Figure  13.  Maximum  value  of  the 
stream  function  versus  Ra/Ract 
from  the  numerical  results  of  this 
report. 


tion  cell  is  smaller  than  that  at  the  onset,  which  explains  the  asymmetrical  behavior  observed  in 
Figure  1 1 . 

Figure  12  shows  the  Nusselt  number  versus  Ra/RaCf  for  each  of  the  same  combinations  of 
boundary  conditions.  A  key  conclusion  from  the  present  results  is  that  simple  relations  for  heat 
transfer  as  a  function  of  Ra/Racr  depend  only  on  the  bottom  thermal  boundary  condition.  For 
an  isothermal  bottom,  Elder’s  equation  (eq  43)  works  well. 

When  the  bottom  is  a  constant  flux  thermal  boundary  condition,  we  propose 


Nu  =  1  +  0.365 


(44) 


as  a  simple  model  of  heat  transfer.  Ribando’s  results  for  the  constant  flux  case  are  also  presented 
in  Figure  12,  for  comparison.  His  results  agree  well  with  eq  44,  although  he  used  an  aspect  ratio 
of  1 .0  for  each  run.  The  actual  aspect  ratio  above  the  onset  of  convection  is  not  known  with  cer¬ 
tainty,  but  we  do  know  that  it  is  smaller  than  that  on  the  onset  of  convection.  Thus„at  least, 
the  geometry  used  by  Ribando  moves  us  in  the  correct  direction  in  approaching  the  optimal  cell 
size.  Similar  results  are  found  with  different  geometries,  so  it  appears  that  the  aspect  ratio  is  of 
little  importance,  as  long  as  the  aspect  ratio  used  is  near  that  at  which  heat  transfer  is  maximized. 

Figures  13  and  14  show  4>max  and  wmax,  respectively,  again  as  functions  of  Ra/Racr.  We 
note  that  again  the  relations  are  dominated  by  the  lower  thermal  boundary  condition.  The  up- 


Ra/ 

'ROcr 

Figure  14.  Maximum  vertical  velocity  versus  Ra/Ract  from 
the  numerical  results  of  this  report. 


per  permeability  boundary  condition  is  more  important  in  determining  the  maximum  velocity,  prob¬ 
ably  because  the  aspect  ratios  change  with  the  upper  boundary  condition.  If  ^max  is  similar,  di/z/dx 
will  decrease  as  the  aspect  ratio  increases.  Thus  it  is  expected  that,  for  a  given  Ra/Racr,  the  max¬ 
imum  velocity  would  be  less  for  an  open  top  layer.  The  values  of  velocity  are  important  in  our 
discussions  of  the  effects  of  convection  on  snow  metamorphism. 

Effects  of  phase  change  on  convection 

Comparing  eq  26  and  29,  we  find  that  the  effect  of  vapor  transport  per  se  arises  from  the  net 
rate  of  phase  change  between  solid  and  fluid.  Henceforth,  we  refer  to  vapor  transport  effects  as 
phase  change  effects,  since  only  with  phase  change  is  there  an  effect  on  the  energy  equation.  As 
shown  in  eq  29,  phase  change  comes  about  as  a  result  of  the  divergence  of  the  vapor  flux.  There 
are  two  contributions  to  the  vapor  flux,  diffusive  and  convective.  Earlier  we  discussed  the  effects 
of  diffusive  vapor  transport  through  snow,  and  found  some  uncertainty  in  the  value  of  De ^  pri¬ 
marily  because  of  uncertainty  over  the  role  of  interparticle  vapor  flow.  Thus,  we  will  vary  Def{ 
to  examine  what  effects  its  magnitude  might  have. 

Earlier  we  defined  the  Rayleigh  number  (eq  5)  as 

po0gATHK 

Ra  =  — - 

UK 

where  k  was  defined  as  km/(pCp)^.  In  our  derivation  of  the  equation  of  motion,  k  was  introduced 
as  a  scaling  factor,  and  could  conceivably  have  several  definitions.  However,  the  energy  equation 
is  greatly  simplified  when  k  is  defined  as  above.  For  a  single-component  fluid  flowing  through  a 
porous  medium,  the  medium’s  thermal  conductivity  is  assumed  to  be  constant  throughout.  For 
snow,  we  have  an  effective  thermal  conductivity  that  strongly  depends  on  temperature,  and  thus  is 
not  suitable  for  use  as  a  scaling  factor.  Thus  we  will  use  km ,  which  is  the  thermal  conductivity  as- 
cribable  to  conduction  through  the  ice  and  air,  and  which  is  essentially  independent  of  temper¬ 
ature.  Note  that  km  =  ke ^  -  kv.  When  the  temperature  of  the  top  surface  is  low  ( T  <  -20°C), 
the  vapor  density  is  very  low  and  thus  transport  of  heat  by  the  vapor  is  negligible,  and  k  used  in 
the  Rayleigh  number  is  truly  descriptive  of  the  heat  transfer  process.  However,  near  the  bottom, 
where  the  temperature  is  higher,  and  vapor  density  higher, «  and  thus  the  Rayleigh  number  do 
not  fully  describe  the  heat  transfer  process. 
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Figure  15.  Ratio  of  heat  transfer  with  phase  change  to  heat  transfer 
without  phase  change  for  two  values  of  the  Rayleigh  number  versus 
Lewis  number,  from  numerical  modeling  of  this  report. 


- With  Phose  Change 


a.  Le  =  0.39. 

Figure  16.  Isotherms  for  convection  with  and  with¬ 
out  phase  change  (Ra  =50;  isothermal,  impermeable 
boundaries;  AR  =  1.0). 

We  can  at  least  qualitatively  analyze  the  effects  of  phase  change,  for  a  given  Rayleigh  number, 
by  weighing  the  relative  contributions  of  the  convective  and  diffusive  phase  change  terms.  If  con¬ 
vective  phase  change  dominates,  then  we  expect  convection  to  be  more  intense  than  analysis  with 
a  single-component  fluid  would  lead  us  to  believe,  and  vice  versa.  We  rewrite  eq  31  as 
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With  Phose  Change 
Without  Phase  Change 
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b.  Le  =  1.0. 

-With  Phose  Change 
-  Without  Phase  Chonge 


c.  Le  =  1.36. 


Figure  16  (cont'd). 


If  the  dimensionless  parameter  in  the  convection  term  is  greater  than  the  parameter  in  the  conduc¬ 
tion  term,  we  expect  convection  to  be  intensified.  This  is  true  if  is  less  than  1.0.  «IDe^ 
is  the  Lewis  number.  If  the  Lewis  number  is  less  than  1 .0,  convection  is  damped,  and  if  it  is  great¬ 
er  than  1.0,  it  is  intensified.  Figure  15  shows  numerical  results  for  various  values  of  the  Lewis 
number.  The  ordinate  is  the  ratio  of  heat  transfer  with  phase  change  to  heat  transfer  without 
phase  change.  The  model  results  closely  follow  the  expected  result  that  near  Le  =  1.0,  phase 
change  has  little  effect  on  heat  transfer.  However,  for  even  small  deviations  from  Le  =  1 .0,  phase 
change  has  significant  effects  on  heat  transfer.  The  range  of  Lewis  numbers  that  the  numerical 
results  cover  gives  approximate  upper  and  lower  limits  for  low  density,  dry  snow. 

Figure  16  shows  the  effects  of  Lewis  number  on  the  temperature  distribution  at  steady  state. 
The  relative  amounts  of  distortion  from  a  linear  temperature  profile  are  more  important  than  the 


exact  positions  of  the  isotherms.  When  the  Lewis  number  is  less  than  1 .0,  the  isotherms  are  less 
distorted  than  the  case  with  no  phase  change,  and  when  the  Lewis  number  is  greater  than  1.0,  the 
isotherms  are  more  distorted.  For  Le  =  1 .0,  the  isotherms  are  of  similar  shape.  The  offset  be¬ 
tween  isotherms  is  a  result  of  the  diffusive  term  changing  the  temperature  profile  in  the  absence 
of  convection. 

Convection  in  sloped  layers 

In  the  Layering  and  Slope  Effects  section,  the  topic  of  convection  in  sloped  layers  was  intro¬ 
duced  and  it  was  noted  that  some  confusion  existed  in  the  literature  over  the  effects  of  slope. 

The  aim  here  is  first  to  clarify  the  role  of  lateral  confinement  on  convection  in  a  sloped  layer, 
and  then  to  quantify  the  effects  of  an  open  top  on  convection  in  a  sloped  layer.  All  discussions 
in  this  section  ignore  the  effects  of  phase  change  and  isolate  the  role  of  slope. 

The  modeling  efforts  here  are  directly  comparable  to  the  experiments  of  Kaneko  et  al.  (1974), 
Their  experiments  were  done  in  a  narrow  slot  of  AR  =  3.0,  which  allowed  only  two-dimensional 
convection.  Their  basic  conclusion  was  that  the  Nusselt  number  could  be  represented  by  the  re¬ 
lation 

Nu  =  0.082(Ra  cos0)0  76  (46) 

tor  10  <  0  <  30  degrees  and  Ra  cos(0)  >  40.  Our  numerical  results  indicate  that  this  is  an  approx¬ 
imate  representation  of  heat  transfer  but  it  obscures  some  of  the  behavior  of  convection  in  an  in¬ 
clined  porous  medium.  Figure  17  shows  the  results  of  Kaneko  et  al.  (1974),  and  the  present  nu¬ 
merical  results  as  a  function  of  Ra  cos(0).  For  a  given  slope,  the  Nusselt  number  increases  in  a 
fashion  similar  to  that  described  by  eq  46,  but  for  a  given  Rayleigh  number,  heat  transfer  either 
increases  or  decreases  with  increasing  slope. 
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Figure  1 7.  Nusselt  number  as  a  function  of  Ra  cos<p  from  this 
report 's  numerical  results  and  the  experimental  results  of  Kan- 
eko  et  al.  ( 1974).  The  envelope  represents  the  extent  of  exper¬ 
imental  scatter  in  the  work  of  Kaneko.  B<uh  our  computations 
and  the  experiments  of  Kaneko  are  dime  for  AR  =  .1.0. 
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Figure  18.  Nusselt  number  from  the  present  numerical  results 
versus  angle  of  inclination.  Computations  done  on  a  grid  of 
aspect  ratio  equal  to  3.0,  with  isothermal,  impermeable  bound¬ 
aries. 


Figure  18  shows  the  Nusselt  number  versus  slope  for  various  values  of  the  Rayleigh  number. 
The  changing  gradient  of  these  curves  can  be  explained  in  terms  of  the  two  competing  forms  of 
convection  that  are  at  work.  We  describe  them  in  terms  of  the  analogous  convection  in  a  fluid 
layer,  the  best  known  of  which  is  the  classic  Benard  convection  in  a  horizontal  layer  with  cold 
top  and  hot  bottom.  Convection  is  multicellular  above  a  critical  value  of  the  Rayleigh  number 
and  increases  in  intensity  as  the  Rayleigh  number  increases.  The  second  form,  or  Rayleigh  con¬ 
vection,  refers  to  flows  driven  by  lateral  temperature  differences.  Convection  then  occurs  for 
any  value  of  the  Rayleigh  number  and  increases  in  intensity  as  the  Rayleigh  number  increases.  It 
is  Rayleigh  convection  that  drives  the  simple  unicellular  flow  described  by  eq  7.  From  eq  7  we 
can  see  that  the  intensity  increases  with  <t>. 

One  of  the  assumptions  leading  to  the  derivation  of  eq  7  is  that  the  temperature  profile  remains 
linear,  or  analogously,  that  heat  transfer  normal  to  the  flow  is  not  affected.  This  is  strictly  true 
only  in  the  limit  as  the  layer  length  becomes  large.  In  a  laterally  confined  layer,  flow  is  curved  by 
end  effects,  and  thus  a  mechanism  for  heat  transfer  between  hot  and  cold  exists.  The  intensity 
of  this  heat  transfer  process  must  increase  as  the  angle  increases  because  the  velocities  increase. 

Weber  (1975)  analytically  treats  the  case  of  convection  in  an  infinite  layer,  and  finds  that  heat 
transfer  is  indeed  a  function  of  Ra  cos0.  In  this  case,  Benard  convection  is  the  only  heat  transfer 
mechanism,  and  thus  the  driving  force  for  Benard  convection  must  decrease  as  the  angle  increases. 

Referring  back  to  Figure  18,  we  see  that  generally  the  Nusselt  number  first  increases  with  in¬ 
creasing  slope,  then  decreases,  and  then  increases  again.  This  happens  because  the  interaction  be¬ 
tween  Rayleigh  and  Benard  convection  results  in  a  changing  cell  pattern.  While  the  heat  transfer 
because  of  Rayleigh  convection  increases  continuously  with  increasing  slope,  the  heat  transfer  be¬ 
cause  of  Benard  convection  decreases  substantially  when  the  number  of  convective  cells  decreases. 
For  a  Rayleigh  number  of  100.0,  the  number  of  cells  decreases  from  three  cells  at  a  slope  of  15° 
to  one  cell  at  a  slope  of  25°.  The  heat  transfer  and  the  Nusselt  number  decrease  accordingly.  At 
the  intermediate  value  of  20°  there  are  still  three  cells  but  the  middle  cell  is  very  weak  and  does 
not  contribute  much  to  the  heat  transfer.  Accordingly,  the  Nusselt  number  is  less  at  20°  than  at 
1 5° .  but  the  big  decrease  occurs  at  25°  because  only  one  cell  is  active.  This  leads  to  the  conclu¬ 
sion  that  while  B'-nard  convection  is  a  more  powerful  heat  transfer  mechanism,  Rayleigh  convec¬ 
tion  also  contributes  to  heat  transfer  when  the  layer  is  confined.  Rayleigh  convection  is  especial¬ 
ly  important  at  Rayleigh  numbers  near  Rarr,  where  Benard  convection  is  weaker. 

If  there  is  no  transition  from  multicellular  to  unicellular  convection,  no  drop  in  heat  transfer 
is  expected  at  slopes  up  to  45°.  Figure  1 9  shows  the  results  of  runs  done  with  an  aspect  ratio  of 
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Figure  19.  Nusselt  number  from  the  present 
numerical  results  versus  angle  of  inclination. 
Computations  done  on  a  grid  of  AR  =1.0, 
with  isothermal,  impermeable  boundaries. 
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Figure  20.  Nusselt  number  versus  Rayleigh  number 
for  two  different  aspect  ratios  of  an  angle  of  inclina¬ 
tion  of  90°.  Our  numerical  results  are  compared  with 
those  of  Chan  etal.  (1970). 


1.0,  where  unicellular  convection  is  forced  on  the  system.  No  transition  in  the  range  of  10° -30° 
is  observed,  but  a  decrease  in  heat  transfer  at  large  (0  >  60°)  angles  is  observed.  This  is  because 
the  driving  force  for  convection  can  be  thought  of  as  the  sum  of  the  Benard  and  Rayleigh  contri¬ 
butions  and  is  thus  proportional  to  cos0  +  sin0.  This  quantity  is  a  maximum  at  45°,  and  we  can 
see  that  at  the  larger  Rayleigh  numbers,  this  is  about  where  the  maximum  occurs.  At  lower  Ray¬ 
leigh  numbers,  Benard  convection  is  not  so  important  in  the  heat  transfer  process,  and  thus  the 
maximum  is  shifted.  As  shown  in  Figure  20,  our  results  agree  well  with  those  in  the  literature  for 
0  =  90°. 

Secondary  shifts  occur  in  Figure  18  at  5°  and  10°  for  Rayleigh  numbers  of  60.0  and  80.0  re¬ 
spectively.  These  small  dips  take  place  when  the  layer  is  in  the  process  of  shifting  cell  forms.  In 
each  case  one  smaller  cell  formed  in  the  middle  of  two  larger  ones.  Although  the  solution  did 
converge,  this  cell  form  would  seem  to  be  unstable,  and  is  probably  an  inefficient  convective  pat¬ 
tern  for  transferring  heat. 

We  now  turn  to  the  important  issue  of  convection  in  a  sloped  layer  of  snow,  where  lateral  con¬ 
finement  is  less  likely  to  be  important.  If  the  snow  layer  is  homogeneous,  but  has  an  impermeable 
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crust  on  its  surface,  convection  at  subcritical  Rayleigh  numbers  is  described  by  eq  7.  If  imperme¬ 
able  layers  are  present  internally  in  the  snow,  eq  7  describes  the  convection  between  the  layers, 
with  the  Rayleigh  number  based  on  the  properties  and  boundary  conditions  of  the  internal  layer. 

Frequently,  the  top  of  the  snow  layer  is  open  to  the  air,  and  the  boundary  condition  is  one  of 
no  flow  in  the  x  direction  (Ribando  1977).  By  following  Combarnous  and  Bories’  (1975)  deriva¬ 
tion  of  eq  7,  we  can  derive  a  similar  expression  for  subcritical  flow  in  a  layer  with  an  open  top. 
Assuming  that  the  layer  is  of  infinite  length,  w  =  0  and  dw/dz  =0, 

^ ~  Tbot  ~  (Tbot  ~  Tt0p)fj  (47) 


where  Ttop  and  Tbot  refer  to  the  temperature  of  the  isothermal  boundaries.  For  an  incompres¬ 
sible  fluid,  the  equation  of  continuity  is  then 


(48) 


The  one-dimensional  form  of  Darcy’s  law  is 


U{z)  =  ~~ufx  '  m  /vK^)U-0(7--ro)]. 


Taking  the  derivative  with  respect  to  z,  we  have 
3 u  K  3 ip  K  ,  .  v  37" 

aT=-^  +  7 **s"*)aF  • 


(49) 


(50) 


This  is  simplified  for  our  case  because  bp/bz  ~  0  and  bT/bz  *  -AT/H,  where  AT  is  Tbot  -  Ttop. 
Integrating  with  respect  to  z,  and  applying  «(//)  =  0,  we  get 


K  AT 

u(z)  = --p0g^(%\n<t>)-rr{z-H).  (51) 

/i  H 

Introducing  the  scaling  factors  x///  and  H  for  velocity  and  distance,  respectively,  we  redefine  u 
and  z  so  that  the  nondimensional  result  is 


m(z)  =  Ra(sin0)  (1-z). 


(52) 


Figure  21  shows  this  analytical  result  along  with  numerical  results  for  a  sloped  layer  of  aspect 
ratio  6.0.  The  agreement  is  good  at  low  angles  of  inclination,  but  deviates  significantly  at  higher 
angles.  This  discrepancy  arises  because  end  effects  are  included  in  the  numerical  model  but  not  in 
the  analytical  model. 

At  higher  Rayleigh  numbers,  there  may  be  multicellular  convection  in  layers  with  an  open  top. 
Given  the  results  of  the  previous  section,  we  expect  that  convection  in  such  a  layer  will  behave  in 
a  manner  similar  to  a  layer  with  a  closed  top.  Thus  we  expect  the  critical  Rayleigh  number  for 
the  onset  of  Benard  convection  in  a  layer  with  permeable  top  and  isothermal  bottom  to  be 


=  =  27 A 

cr  COS 0  COS0 


(53) 
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Figure  21.  Velocity  profiles  for  convection  in  a  sloped  layer  with  a 
permeable  top.  Analytical  (eq  51)  and  numerical  results  for  K  = 

2.4x  10~4  cm2,  AT  =  20° C.  Analytical  results  are  for  an  infinitely 
long  layer,  while  numerical  results  are  for  AR  =  6.0. 

We  also  expect  heat  transfer  results  to  follow  the  same  patterns  as  in  a  closed  layer,  with  Benard 
convection  dominating  in  infinite  layers  and  both  Rayleigh  and  Benard  convection  influencing 
heat  transfer  in  confined  layers. 

A  numerical  model  has  been  developed  and  utilized  to  describe  thermal  convection  in  snow. 

In  particular,  the  influences  of  constant  flux  bottom  and  permeable  top  boundary  conditions,  of 
phase  change  and  of  inclination  of  the  snow  cover  have  been  described  and  quantified.  Each  is 
found  to  have  a  significant  effect  on  both  the  intensity  and  frequency  of  convection.  The  im¬ 
portant  topic  of  layering  has  not  been  discussed  here,  but  has  been  dealt  with  in  the  literature  and 
is  well  understood.  In  the  following  sections,  we  discuss  experiments  performed  in  an  effort  to 
further  test  our  numerical  model,  and  then  apply  the  findings  of  this  section  in  discussing  the  im¬ 
portance  of  convection  in  snow  covers. 


EXPERIMENTS 

Introduction 

We  now  experiment  with  convection  through  a  horizontal  porous  layer  that  is  confined,  with 
a  warm  bottom  and  cool  top.  Similar  experiments  have  frequently  been  done  by  others  working 
with  porous  media,  but  some  significant  differences  exist.  Here  the  lower  thermal  boundary  con¬ 
dition  is  a  constant  flux,  rather  than  a  constant  temperature  boundary  as  has  been  used  previously. 
In  addition,  when  snow  is  used,  phase  change  effects  are  an  important  component  of  the  heat 
transfer. 

We  will  describe  the  apparatus  used,  the  data  acquisition  equipment  and  the  data  analysis. 

The  experimental  results  are  presented  in  the  next  section  and  we  compare  these  results  with  rel¬ 
evant  numerical  and  theoretical  results. 

Experimental  apparatus 

Experiments  were  done  using  glass  beads  with  air  as  the  fluid,  glass  beads  with  water  and  dry 
snow  with  air.  The  apparatus  used  in  the  glass  bead  experiments  was  substantially  different  from 
that  used  in  the  snow  experiments,  but  the  data  collection  system  was  essentially  the  same. 


Figure  22.  Apparatus  far  glass  bead  experiments. 


a.  Snow  samples  separated  by  insulation. 

Figure  23.  Apparatus  for  snow  experiments. 

Experiments  were  first  performed  on  glass  beads  to  gain  experience  in  doing  experiments  on 
natural  convection  before  moving  to  the  more  difficult  and  unknown  problems  of  convection  in 
snow.  Figure  22  shows  the  apparatus  used  for  the  glass  bead  experiments.  At  the  top  boundary 
cold  air  was  blown  over  an  aluminum  plate;  the  air  temperature  was  kept  constant  by  using  a  heat¬ 
er,  connected  to  a  temperature  controller.  At  the  bottom  boundary,  electric  heat  pads  were  used 
to  supply  the  constant  heat  flux.  The  flux  rate  was  controlled  by  regulating  the  voltage  input  us¬ 
ing  a  variable  power  supply.  The  glass  bead  layer  was  approximately  40  by  40  cm  by  15  cm  deep. 
When  water  was  the  saturating  fluid,  a  plastic  liner  was  used  to  retain  it. 

Figure  23a  shows  a  three-dimensional  view  of  the  snow  experiment,  while  Figures  23b  and  c 
are  schematics  of  the  experiment.  Because  of  the  low  thermal  conductivity  of  snow,  heat  losses 
out  the  sides  would  have  been  unacceptably  large  if  the  same  apparatus  was  used  for  both  glass 
beads  and  snow.  A  guarded  technique  similar  to  that  of  Buretta  and  Berman  (1976)  was  used  with 
a  tall,  narrow  snow  sample  (50  by  40  cm  by  10  cm  high)  surrounded  by  polystyrene  insulation. 

The  large  sides  (i.e.,  front  and  back)  had  identical  snow  samples  on  the  other  side  of  the  insulation. 
By  heating  the  bottoms  and  cooling  the  tops,  we  could  maintain  similar  temperature  profiles  in¬ 
side  and  outside  the  insulation,  thus  reducing  the  driving  forces  for  heat  loss  from  the  central  snow 
sample.  The  narrow  ends  had  a  loose  insulation  composed  of  styrofoam  packing  materials  outside 
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b.  Side  view. 


c.  Front  view. 


Figure  23  f cant ’d).  Apparatus  for  snow  experiments. 


the  polystyrene  slabs.  This  loose  insulation  was 
also  cooled  on  the  top  and  heated  on  the  bottom. 

Each  heater  was  wired  to  a  separate  voltage  source, 
so  that  temperature  profiles  in  each  snow  sample 
could  be  matched  as  closely  as  possible.  However,  Temperoture 
only  the  central  snow  sample  was  used  for  data 
collection. 

The  top  boundary  for  the  snow  experiments 
was  an  aluminum  cooling  plate  through  which 
cold  glycol  flowed  (Fig.  24).  In  essence,  the  cool¬ 
ing  plate  was  a  rectangular  box,  with  internal  di¬ 
mensions  of  about  50  by  50  cm  by  2  cm  deep. 

The  glycol  was  cooled  by  a  copper  coil  through 
which  the  refrigerant  flowed.  The  glycol  was 
kept  at  a  constant  temperature  by  using  a  small  Figure  24.  Glycol  supply  system  for  the 
heater  with  an  on-off  temperature  controller.  snow  experiments. 
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Figure  25.  Temperature  measurement  apparatus 
and  thermocouple  placement. 


Control  of  the  glycol  temperature  was  good  enough  to  eliminate  cyclic  temperature  variations  at 
the  top  snow  surface. 

The  temperature  measurement  system  (Fig.  25)  consisted  of  33  copper-con stantan  (Type  T) 
thermocouples,  wired  to  a  36-position  switch  that  fed  a  Fluke  digital  thermometer  (accuracy 
±0.1  °C).  Thermocouples  were  distributed  in  the  central  snow  sample  as  illustrated  in  Figure  25. 
Their  accuracy  was  verified  using  an  ice-water  bath  between  each  experiment.  In  each  check,  the 
thermocouples  read  either  -0.1°  or  -0.2°C.  Probably  a  slight  offset  existed  in  the  Fluke  but  the 
offset  was  not  significant  since  we  were  mostly  interested  in  temperature  differences.  The  uncer¬ 
tainty  in  the  thermocouple  measurements  was  taken  as  ±0.1  °C,  which  is  negligible. 

From  the  experimental  data  we  calculated  the  net  vertical  heat  flux  through  and  effective  ther¬ 
mal  conductivity  of  each  sample.  A  plot  of  effective  thermal  conductivity  versus  net  heat  flux 
has  the  same  physical  meaning  as  a  plot  of  Nusselt  number  versus  Rayleigh  number  when  the  Ray¬ 
leigh  number  is  based  on  a  flux  boundary  condition.  The  dimensional  form  is  used  for  clarity. 

In  these  experiments  the  temperature  profiles  and  the  heat  input  at  the  bottom  boundary  were 
measured  directly.  If  no  heat  was  lost  through  the  boundaries  of  the  experiment,  then  the  net 
vertical  heat  flux  was  simply  the  heat  input  per  unit  area.  However,  heat  did  flow  through  the 
sidewalls.  These  lateral  heat  flows  must  be  accounted  for  in  the  calculation  of  the  net  vertical 
heat  flux  but  the  accuracy  of  the  computation  is  limited  by  the  knowledge  we  have  of  the  temper¬ 
ature  distributions  inside  and  outside  the  sidewalls. 

Lateral  heat  flows  were  calculated  using  the  temperature  difference  across  the  midpoint  of  the 
insulation.  If  the  vertical  temperature  profiles  inside  and  outside  the  sidewalls  are  both  linear, 
then  the  temperature  difference  across  the  sidewall  also  varies  linearly  with  height  so  that  the  aver¬ 
age  temperature  difference  is  that  at  the  midpoint.  Unfortunately,  we  have  no  means  to  deduce 
how  close  the  profiles  were  to  linear  and  thus  the  maximum  value  of  the  uncertainty  in  the  heat 
flux  calculation  was  taken  to  be  the  total  calculated  heat  flux  through  the  sidewalls. 
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RECOMMENDATIONS 


In  this  work,  the  means  by  which  a  fairly  complete  understanding  of  convection  in  snow  might 
be  achieved  are  developed.  However,  definitive  conclusions  cannot  be  made  because  we  lack  a 
good  quantitative  understanding  of  many  fundamental  processes  in  and  properties  of  snow. 

f  oremost  among  these  ups  is  a  lack  of  knowledge  of  just  how  much  vapor  is  transferred  through 
snow  by  diffusion.  In  essence,  we  do  not  know  the  macroscopic  (effective)  diffusivity  of  water 
vapor  through  snow.  Simple  experiments  done  20  or  30  years  ago  constitute  our  entire  knowledge. 
These  results  supply  hints,  but  in  many  ways  seem  inconsistent  with  our  present  understanding 
of  vapor  transfer  through  snow. 

Related  to  our  lack  of  understanding  of  mass  transfer  through  snow,  but  with  further  compli¬ 
cations  of  its  own,  is  the  subject  of  heat  transfer.  At  present,  our  knowledge  of  thermal  conduc¬ 
tivity  is  limited  to  widely  scattered  data  on  effective  thermal  conductivity,  which  lump  together 
conduction  through  ice,  conduction  through  air  and  latent  heat  transfer  due  to  vapor  diffusion. 

Our  understanding  of  the  relative  contribution  of  each  component  is  purely  speculative.  A  thor¬ 
ough  understanding  ot  how  heat  is  transferred  through  snow  is  important  not  only  in  macroscop¬ 
ic  treatments  of  snow  like  this,  but  also  in  microscopic  descriptions  of  the  snow  metaniorphism 
process. 

Some  understanding  of  the  macroscopic  diffusivity  of  snow  could  be  achieved  by  measuring 
effective  thermal  conductivity  as  a  function  of  temperature  only.  In  addition,  a  knowledge  of 
the  tortuosity  of  the  ice  lattice  can  be  gained  by  measuring  the  electrical  resistance  of  snow.  By 
coordinating  the  two  efforts,  it  is  felt  that  a  better  understanding  of  heat  transfer  through  snow 
is  possible.  Calculation  of  diffusion  coefficients  might  also  be  possible  with  some  recent  models 
of  local  vapor  transport. 

Our  knowledge  of  snow  properties  in  field  situations  could  also  be  improved.  In  particular, 
it  would  be  interesting  to  see  thermal  conductivity  and  permeability  measurements  done  at  the 
same  field  sites.  This  could  give  us  a  much  better  idea  of  the  ways  the  two  are  related,  and  allow 
a  more  definite  statement  concerning  the  frequency  with  which  convection  takes  place. 

As  was  indicated  earlier,  these  experiments  on  convection  in  snow  should  be  extended  to  test 
the  initial  results.  It  seems  important  to  measure  the  permeability  of  the  sample,  so  that  more 
definite  conclusions  can  be  made  regarding  the  agreement  between  experiment  and  theory.  In 
addition,  sensitive  experiments  are  needed  to  confirm  the  validity  of  theories  on  the  role  of 
phase  change  on  convection.  The  experimental  difficulties  of  doing  very  accurate  heat  transfer 
experiments  on  snow  are  apparent  and  such  experiments  will  require  much  time  and  patience. 

finally,  experiments  aimed  at  understanding  the  role  of  convection  in  snow  metamorphism 
are  needed,  along  with  better  theoretical  treatments.  A  fairly  simple  experiment  would  be  to  tilt 
a  sample  of  snow,  so  that  convection  is  forced  on  the  system.  In  this  manner,  both  convection 
perpendicular  to  the  temperature  gradient  and  convection  parallel  to  the  gradient  could  be 
studied. 
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SUMMARY 


The  purpose  of  this  work  was  to  investigate  thermally  driven  flows  through  snow.  Toward 
this  end,  a  numerical  model  was  constructed  and  experiments  were  performed  on  snow  and  glass 
beads. 

The  starting  point  of  the  investigation  was  a  theory  developed  to  describe  flows  through  po¬ 
rous  media  in  general.  Extensions  of  the  theory  were  made  using  a  numerical  model  that  accounts 
for  boundary  conditions  of  particular  interest  for  the  study  of  snow.  Simple  relations  were  de¬ 
veloped,  based  on  the  results  of  the  model,  which  can  be  used  to  predict  heat  transfer  based  on 
the  choice  of  the  lower  thermal  boundary  condition  and  the  specification  of  the  parameter  Ra/ 
Racf..  The  relations  are  valid  for  either  permeable  or  impermeable  upper  boundaries.  For  the 
case  of  an  isothermal  lower  boundary,  the  results  fit  existing  experimental  data  as  well  or  better 
than  analytical  results.  For  a  constant  tlux  lower  boundary,  the  results  of  our  experiments  agree 
well  with  the  present  numerical  calculations. 

Latent  heat  terms  accounting  for  phase  change  between  water  vapor  and  ice  are  included  in  a 
unique  energy  equation  used  in  the  numerical  model.  These  terms  were  found  to  be  significant, 
especially  at  the  low  Rayleigh  numbers  of  interest  for  snow.  The  phase  change  term  can  intensi¬ 
ty  or  weaken  convection,  depending  on  the  ratio  of  the  macroscopic  diffusivity  to  the  thermal 
conductivity  ot  the  snow.  Phase  change  terms  also  have  a  small  effect  on  the  conditions  for  the 
onset  of  convection. 

Convection  in  sloped  layers  was  also  examined  numerically.  It  was  found  that  misunderstand¬ 
ings  existed  in  previous  works  dealing  with  convection  in  a  sloped,  porous  layer.  Here  it  is  found 
that  the  intensity  of  convection  is  not  just  a  function  of  the  parameter  Ra  cos 0  but  that  lateral  con¬ 
finement  plays  a  very  significant  role.  A  physical  explanation  is  offered  that  consistently  explains 
the  observed  results  in  terms  of  a  transition  from  Benard  type  convection  to  Rayleigh  type  convec¬ 
tion. 

Analytical  solutions  were  developed  to  describe  the  velocity  field  for  convection  in  a  sloped 
layer  with  a  permeable  top  at  low  Rayleigh  numbers.  The  analytical  solution  agrees  well  with  the 
numerical  model.  The  importance  of  layering  on  convection  in  sloped  layers  at  low  Rayleigh  num¬ 
bers  is  qualitatively  explained. 

Experiments  indicate  that  there  is  natural  convection  through  snow.  Because  of  uncertainties 
in  the  permeability,  it  was  necessary  to  calibrate  the  modeling  by  inferring  the  critical  heat  flux 
value  Horn  the  observed  onset  of  convection.  Except  for  this  one  problem,  agreement  between 
theory  and  experiment  was  good.  The  permeability  implied  by  the  critical  heat  flux  we  selected 
is  different  from  that  which  we  would  expect  based  on  the  few  reported  values  in  the  literature. 

This  points  to  an  obvious  need  for  further  experimental  work  on  the  physical  properties  of  snow, 
especially  permeability. 

The  existing  theory  for  the  onset  of  convection  in  porous  media  was  applied  to  two  different 
field  settings.  We  found  that  convection  is  unlikely  in  the  alpine  regions  of  the  North  American 
Rockies,  but  is  likely  to  be  very  common  in  subarctic  regions  such  as  interior  Alaska.  Simple,  uni¬ 
cellular  convection  in  sloped  snowpacks  is  expected  in  both  regions. 

Finally,  simple  theories  to  describe  the  effects  of  convection  on  snow  metamorphism  were  in¬ 
troduced  and  applied.  For  velocities  calculated  from  the  modeling  results,  we  found  that  convec¬ 
tion  should  have  little  effect  on  snow  crystal  growth.  However,  convection  was  found  to  have  a 
substantial  effect  on  the  transfer  of  vapor  through  the  snow  cover.  It  seems  possible  that  thermal 
convection  may  not  substantially  increase  the  rate  of  metamorphism  but  might  trigger  the  growth 
of  the  more  spectacular  hollow  crystals  with  scrolls  and  striations. 
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Intuitively,  it  would  seem  that  almost  any  convection  would  have  an  important  impact  on  the 
transfer  of  mass,  given  the  slow  rate  at  which  vapor  is  transferred  by  diffusion.  A  simple  method 
of  comparing  macroscopic  convective  and  diffusive  fluxes  is  by  examining  the  ratio 


n  dPv dT 
ueff - 

a  T  dz 

For  common  values  of  dT/dz  of  0.2°C/cm  and  v  of  0.04  cm/s,  this  ratio  is  about  10,  indicating 
that  convection  is  very  important  in  the  overall  mass  transfer  process.  This  analysis  says  nothing 
about  what  happens  to  individual  grains,  only  about  how  much  mass  is  transferred  through  the 
snow. 

The  previous  discussions  indicate  that  while  convection  is  important  in  moving  vapor,  it  only 
moves  vapor  past  particles  and  has  little  impact  on  crystal  growth.  Thus,  even  when  convection 
occurs,  snow  metamorphism  is  controlled  by  the  rate  of  local  mass  transfer.  However,  it  seems 
unreasonable  to  conclude  that  an  order  of  magnitude  increase  in  macroscopic  mass  flux  leads  to 
no  increase  in  local  mass  flux,  so  perhaps  our  means  of  quantifying  mass  transfer  to  a  particle  are 
inadequate. 

The  above  methods  for  describing  the  rates  of  mass  transfer  to  a  particle  under  convective-dif¬ 
fusion  conditions  assume  that  the  flow  near  the  particle  is  entirely  forced  convection,  that  is,  that 
the  transfer  processes  near  the  particle  have  little  effect  on  the  flow  field.  Donaghey  (1980)  es¬ 
tablished  the  following  criteria  for  the  transition  from  forced  to  mixed  free  convection: 


Gr/Re2  <  0.3 


forced  convection 


0.3  <  Gr/Re2  <  16  mixed  convection 


16  <  Gr/Re2 


free  convection 


where  Gr  is  the  Grashof  number,  defined  as  p2PgAT L3 Ip.2 .  By  defining  the  characteristic  length 
/.  as  one  particle  diameter  and  defining  the  temperature  difference  as  this  length  times  an  average 
temperature  gradient,  we  find  that  the  Grashof  number  is  about  0.055.  For  a  velocity  of  0.04  cm/s, 
the  Reynolds  number  is  about  0.03,  and  thus  the  ratio  Gr/Re  is  about  6.  This  indicates  that  even 
over  the  limited  range  of  one  crystal  diameter,  buoyancy  forces  are  important  in  adequately  describ¬ 
ing  the  flow  field. 

To  fully  analyze  the  effects  of  convection  on  snow  metamorphism,  it  would  seem  to  be  necessary 
to  solve  the  fully  coupled  mass,  energy  and  momentum  equations  around  a  particle.  Because  of  the 
complex  geometries  and  mathematics,  this  is  a  formidable  task,  which  can  be  addressed  only  by 
numerical  methods. 

Finally,  we  note  the  experimental  work  of  Keller  and  Hallet  (1982)  on  ice  crystal  growth  in  a 
convective  environment.  Although  the  supersaturations  and  velocities  used  were  several  orders  of 
magnitude  greater  than  what  we  expect  in  snow,  their  results  may  be  applicable.  A  key  conclu¬ 
sion  of  theirs  was  that  even  in  the  lower  ranges  of  velocities  investigated,  the  introduction  of  any 
flow  could  change  the  crystal  form  from  a  solid  type  (plates,  columns)  to  a  skeletal  type  (den¬ 
drites,  needles).  This  transition  took  place  even  when  the  change  in  the  growth  rate  was  small. 

They  attributed  the  skeletal  forms  to  enhanced  vapor  density  gradients  along  the  crystal  face, 
which  leads  to  more  rapid  growth  near  corners.  This  suggests  that  the  onset  of  convection  in  some 
way  may  coincide  with  the  onset  of  the  growth  of  hollow  crystals. 


is  n 
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Applications  to  snow  metamorphism 

Snow  metamorphism  is  a  process  of  crystal  growth,  the  complex  interaction  of  heat  removal 
from  and  mass  transfer  to  a  crystal  with  the  kinetics  of  phase  change  at  the  crystal  surface.  In 
this  section,  we  concentrate  on  the  mass  transfer  process,  that  is,  on  the  means  by  which  vapor 
is  transported  from  vapor  source  to  sink,  introducing  various  theoretical  means  of  assessing  the 
importance  of  convection  to  the  mass  transfer  process.  This  section  is  really  more  of  a  literature 
review  than  new  work;  its  intent  is  to  point  the  way  for  future  workers. 

Colbeck  (1983)  assesses  the  effects  of  convection  in  terms  of  a  ventilation  factor,  the  amount 
by  which  the  mass  growth  rate  is  increased  because  of  flow.  This  approach  was  developed  by 
cloud  physicists  as  a  means  of  assessing  the  influence  of  fall  velocity  on  crystal  growth  rates  (see 
Prupacher  and  Klett  [1978]  for  a  thorough  review).  Although  results  vary  somewhat,  depending 
on  which  correlation  for  the  ventilation  factor  is  used,  in  general  it  is  safe  to  say  that  only  at  ve¬ 
locities  of  about  1  cm/s  or  greater  will  the  influence  be  substantial. 

From  Figure  14  we  can  see  that  a  typical  dimensionless  maximum  velocity  for  Rayleigh  num¬ 
bers  just  above  the  critical  value  is  about  4.0.  The  scaling  factor  is  k/H,  or  km  /( pCp  )^N,  where 
typical  values  are  km  of  2.5x1  O'4  cal/cm  s  °C,  of  1 ,36x  1 0"’  g/cm3 ,  Cp  of  0.24  cal/g  °C  and 
//  of  100  cm.  Then  the  dimensional  velocity  is  w  =  0.04  cm/s.  This  is  approximately  the  same 
velocity  that  results  from  solving  eq  52  for  Ra  =  10.0  and  an  angle  of  30°,  and  thus  is  a  good 
choice  as  a  typical  velocity  for  use  in  estimating  the  effects  of  either  Benard  or  Rayleigh  convec¬ 
tion. 

Obviously,  a  velocity  of  0.04  cm/s  is  well  below  that  which  the  ventilation  factor  approach 
suggests  is  important  in  mass  transfer  to  a  particle.  However,  it  is  not  certain  that  the  ventilation 
factor  approach  is  a  valid  one.  The  key  parameter  in  describing  flow  around  a  particle  is  the  Rey¬ 
nolds  number,  defined  as 


Uf 


(57) 


where  v  is  the  undisturbed  velocity.  Reynolds  numbers  for  flow  around  falling  snow  particles  are 
frequently  in  the  range  1 .0-10.0,  large  enough  so  that  a  laminar  boundary  layer  develops.  The 
thickness  of  this  boundary  layer  (which  decreases  with  increasing  Reynolds  number)  controls  the 
rate  of  mass  transfer.  The  boundary  layer  is  best  thought  of  as  the  distance  over  which  mass  must 
diffuse  to  get  from  the  free  stream  concentration  to  the  surface  concentration.  Thus,  as  the 
boundary  layer  thickness  decreases,  the  rate  of  mass  transfer  increases. 

For  thermal  convection  through  snow,  velocities  are  low  enough  so  that  the  flow  may  be  de¬ 
scribed  as  creeping  flow,  in  which  viscous  terms  are  important  everywhere.  The  boundary  layer 
has  no  meaning  in  creeping  flow  situations  (in  effect,  the  boundary  layer  thickness  is  infinite),  and 
thus  the  boundary  layer  solutions  should  tell  little  about  mass  transfer  during  convection  through 
snow. 

The  effects  of  creeping  flow  on  mass  transfer  to  a  single  sphere  have  been  solved  by  Acrivos 
and  Taylor  (1962).  Their  results  indicate  that 

Sh  =  1 +~  +  ^  In  Pe  +0.03404  Pe3 +0(Pe3)  (58) 


for  Pe  <  1 .0,  where  Pe  is  the  Peclet  number,  defined  as  vd/D,  where  d  is  the  particle  diameter  and 
D  the  diffusion  coefficient.  In  a  rough  sense,  the  Peclet  number  is  the  ratio  of  convective  to  dit- 
fusive  contributions.  The  Sherwood  number  Sh  is  similar  to  a  ventilation  coefficient.  It  describes 
the  actual  increase  in  mass  transfer  when  flow  is  added  to  the  diffusion.  For  the  velocities  described 
above,  however,  Peclet  numbers  are  small  and  this  solution  predicts  that  flow  has  little  impact  on 
mass  transfer. 
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isothermal.  It  is  not  known  at  this  point  whether  the  bottom  boundary  condition  is  a  constant 
temperature  or  a  constant  flux  condition,  and  thus  both  will  be  considered.  Applying  the  results 
from  Table  1 ,  we  find  the  critical  temperature  difference 

i30.6°C  constant  flux  bottom 

46.9°C  isothermal  bottom 

While  temperature  differences  of  30°C  may  occur  in  alpine  areas,  they  are  certainly  not  common, 
and  thus  this  analysis  would  indicate  that  Benard  type  convection  is  probably  not  important  in 
alpine  areas.  In  this  analysis  we  have  chosen  values  of  the  thermal  conductivity  and  permeability 
based  on  our  best  estimates  of  field  conditions. 

In  subarctic  regions,  such  as  interior  Alaska,  convection  is  more  likely.  Trabant  and  Benson 
( 1972)  report  that  snow  densities  are  frequently  around  0.2  g/cm3 .  and  that  by  late  winter  the 
snow  cover  can  be  composed  primarily  of  depth  hoar  crystals.  These  may  be  as  large  as  1  cm. 
Typical  snow  covers  are  50-80  cm  deep,  and  are  commonly  subjected  to  temperature  differences 
of  30°C.  Although  we  have  no  permeability  or  thermal  conductivity  data  for  this  region,  it  is 
likely  that  the  ratio  of  permeability  to  thermal  conductivity  is  greater  than  for  alpine  snow  covers 
by  at  least  an  order  of  magnitude.  This  result  indicates  that  convection  may  be  common  in  re¬ 
gions,  such  as  interior  Alaska,  where  large  temperature  gradients  are  imposed  upon  the  snowpack. 
We  take  these  conclusions  as  tentative  pending  the  acquisition  of  better  data  on  permeability  and 
conductivity. 

In  this  discussion,  we  have  not  accounted  for  the  effects  of  phase  change,  layering  or  slope. 
Phase  change  can  raise  or  lower  the  critical  temperature  difference  necessary  for  the  onset  of  Be¬ 
nard  convection.  For  a  conductivity  of  2. Ox  10"4  cal/cm  s  °C,  the  Lewis  number  is  1.0  when  the 
macroscopic  diffusivity  is  0.61  cm2  /s  and  phase  change  has  no  effect  on  the  onset  of  Benard  con¬ 
vection.  If  the  diffusivity  is  higher,  as  Yen  (1963)  and  Yosida  (1955)  have  observed,  the  onset 
of  convection  is  at  a  higher  temperature  difference  than  we  earlier  calculated  and  vice  versa  at 
lower  diffusivities.  For  a  diffusivity  of  0.22  cm2/s,  that  of  water  vapor  through  air,  and  probably 
a  lower  bound,  the  effect  is  to  lower  the  critical  temperature  difference  about  10%.  If  the  diffu¬ 
sivities  found  by  Yen  and  Yosida  are  correct,  the  effect  of  phase  change  on  the  onset  of  convec¬ 
tion  is  to  increase  the  critical  temperature  difference  by  less  than  10%. 

The  effects  of  layering  are  harder  to  quantify  but  are  probably  not  negligible.  However,  in  the 
theoretical  analysis,  we  used  properties  of  snow  that  are  conducive  to  the  onset  of  convection. 
Thus  any  layering  effect  would  have  a  negative  impact  on  convection.  The  exception  to  this  is 
when  there  are  significant  depth  hoar  layers,  in  which  case  the  permeability  of  individual  layers 
can  be  substantially  greater  than  the  permeability  used  above.  If  the  only  layering  is  that  ascrib- 
able  to  the  growth  of  depth  hoar  layers,  then  convection  may  be  more  likely  than  we  have  indi¬ 
cated. 

The  onset  of  Benard  convection  in  sloped  layers  is  given  by  (see  the  Convection  is  Sloped  lay¬ 
ers  section) 


Racr 


COS  0 


where  Rafr  n  is  the  critical  Rayleigh  number  for  a  horizontal  layer.  Thus  Benard  convection  is 
less  likely  in  a  sloped  layer  than  in  a  horizontal  layer.  In  alpine  areas  it  is  very  unlikely;  while  in 
interior  Alaska  it  is  probably  still  common,  since  even  for  slopes  of  30°  the  increase  in  the  critical 
temperature  difference  is  only  about  16%.  However,  it  is  important  to  remember  that  Rayleigh 
convection,  a  single  cell  pattern  flowing  parallel  to  the  slope,  is  important  in  all  regions. 
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APPLICATIONS  AND  CONCLUSIONS 


In  this  section,  the  aim  is  to  use  what  we  have  learned  to  determine  the  consequences  of  con¬ 
vection  in  snow.  The  treatment  is  an  introduction  to  what  may  be  done  with  our  present  under¬ 
standing  of  convection  through  snow. 

We  will  focus  on  two  topics.  First,  the  theoretical  treatment  of  the  onset  of  convection  is  ap¬ 
plied  to  natural  snow  covers  in  an  attempt  to  determine  just  how  common  convection  is  under 
existing  conditions.  Second,  some  approaches  to  modeling  snow  metamorphism  under  convec¬ 
tive  conditions  are  introduced. 


Onset  of  Benard  convection  in  seasonal  snow  covers 

Previously,  we  developed  the  theoretical  basis  for  predicting  the  onset  of  Benard  convection, 
based  on  the  value  of  the  Rayleigh  number  and  the  boundary  conditions.  This  development  al¬ 
lows  the  prediction  of  the  onset  for  either  horizontal  or  sloped  layers.  Rayleigh  convection  will 
occur  for  any  inclined  snow  cover  with  an  applied  gradient.  We  calculated  the  critical  conditions 
for  the  onset  of  convection,  and  then  indicated  how  common  Benard  convection  might  be  under 
naturally  occurring  conditions. 

The  Rayleigh  number  is  defined  in  eq  5.  Properties  of  the  snow  cover  appear  as  the  snowpack 
depth  //,  the  temperature  difference  based  on  boundary  conditions  AT,  the  intrinsic  permeabili¬ 
ty  K  and  the  thermal  conductivity  of  the  porous  medium  km .  The  remaining  parameters  are  tem¬ 
perature-dependent  fluid  properties  or  constants.  For  air,  we  use  properties  at  typical  conditions 
of -15°C  and  10s  Pa  (=^l  atm)  from  Table  5.  The  Rayleigh  number  is  then 


Ra  =  0.0105 


ATHK 


In  the  Effects  of  Phase  Change  on  Convection  section,  we  showed  that  the  correct  thermal 
conductivity  to  use  was  km ,  or  that  based  on  conduction  only,  i.e.,  not  an  effective  conductivity 
that  includes  vapor  diffusion.  We  ignored  both  diffusive  and  convective  phase  change  terms  that 
gave  the  same  results  as  incorporating  both  components  when  the  Lewis  number  is  1 .0,  a  reason¬ 
able  approximation  for  low  density,  dry  snow.  Later  in  this  section  we  will  further  discuss  effects 
of  phase  change  on  the  onset  of  convection. 

Unfortunately,  no  data  are  available  that  clearly  delineate  the  contributions  of  conduction  and 
diffusion  to  the  heat  transfer  process.  In  the  Heat  Transfer  section,  we  estimated  the  contribution 
of  vapor  diffusion  to  the  effective  thermal  conductivity  as  about  2.5x  KT4  cal/cm  s  °C  for  a  value 
of  the  diffusion  coefficient  (proposed  by  Yosida  1955)  of  0.85  cm2/s.  In  alpine  areas  where  the 
growth  of  depth  hoar  is  common,  a  usual  snowpack  density  is  about  0.25  g/cm3  (Giddings  and 
LaChapelle  1962).  Effective  thermal  conductivity  data  compiled  by  Mellor  (1977)  show  that  for 
this  density,  effective  thermal  conductivities  range  from  2.5x  lO'4  to  6.0x  10"4  cal/cm  s  °C.  The 
pure  conduction  conductivity  is  then  probably  around  2.5x  lO"4  cal/cm  s  °C. 

For  an  arbitrary  grain  size  of  1  mm,  and  a  density  of  0.25  g/cm3 ,  Shimizu  (1970)  predicts  a 
permeability  of  1 . 1  x  1  O'4  cm.  For  a  snow  cover  1 00  cm  deep,  the  critical  temperature  difference 
necessary  for  the  onset  of  convection  is 

ATcr  =  1 .73  Racr  °C  . 

The  critical  Rayleigh  number  is  a  function  only  of  boundary  conditions  (ignoring  phase  change 
effects).  The  soil/snow  interface  is  relatively  impermeable  and  the  top  boundary  is  normally  per¬ 
meable,  although  a  wind  or  rain  crust  can  reduce  its  permeability.  The  top  boundary  is  in  con¬ 
tact  with  air  that  is  very  nearly  at  a  constant  temperature,  thus  the  top  snow  boundary  is  nearly 
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a.  Lateral  profiles. 
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b.  Vertical  profiles.  Solid  line  indicates  the 
conduction  solution. 


Figure  30.  Temperature  profiles  for  experiments  with  snow  for  two  runs 
where  convection  occurred  and  one  where  it  did  not. 


Figure  31.  Effective  thermal  conductivity  versus  net  heat  flux. 
Solid  line  represents  results  of  numerical  model  without  phase 
change.  Dashed  line  represents  numerical  results  including 
phase  change  with  Deff  equal  to  0.05  cm*  /s. 


vection  does  occur,  with  a  value  of  critical  heat  llux  of  around  l.Hx  U)4  cal  cm’  s.  Calculation 
of  the  critical  Rayleigh  number  for  the  onset  of  convection  is  difficult  for  snow  ,  primarily  because 
of  uncertainty  in  the  permeability.  Thermal  conductivity  was  measured  in  our  experiments  and 
thus  poses  no  real  problem. 

A  survey  of  the  snow  at  the  start  of  the  experiment  revealed  an  average  grain  si/e  of  about  1 .5 
mm.  The  snow  used  was  from  the  bottom  5  cm  of  the  snow  cover,  beneath  a  thick  crust.  The 
large  grains  were  formed  by  wet  snow  metamorphism  followed  by  dry  metamorphism  because  of 
strong  temperature  gradients.  The  original  snow  layer  was  very  permeable,  as  the  grain  clusters 
formed  by  the  wet  snow  metamorphism  were  still  evident,  although  weakly  bonded  together.  Ob¬ 
servations  of  the  snow  at  the  end  of  the  experiment  indicated  that  there  was  no  significant  grain 
growth  during  the  experiment. 

The  sample  was  prepared  bv  disaggregating  the  grain  clusters  and  carefully  filling  the  box,  so 
as  to  maintain  as  low  a  density  as  possible.  Although  care  was  taken,  the  snow  in  the  box  ap¬ 
peared  much  less  permeable  than  the  layer  outside.  The  density  was  0,25  g/cm3 . 

There  is  much  uncertainty  over  the  intrinsic  permeability  of  snow.  Shimizu  (1970)  proposed 
the  formula 

K  -  0.077c/ :  exp(-7.8  pjpw )  (56) 

for  snow  of  mean  diameter  less  than  1  mm.  Shimizu  measured  the  permeability  of  coarse-grained 
snow  at  densities  of  0.4  g/cnt3  or  higher,  and  found  a  range  of  permeabilities  of  about  2.6x  10's 
to  5.1  x  10"*  cm2 .  Bader  (1939)  measured  permeabilities  for  coarse-grained  snow  and  depth  hoar 
over  the  density  range  of  our  interest.  He  found  0.85x  Iff4  <  K  <  1 ,36x  10-4  cm2  for  coarse¬ 
grained  snow,  and  2. Ox  10^  <  K  <  2.5x  ICf4  cm2  for  depth  hoar.  Shimizu,  however,  believed 
that  Bader's  experimental  technique  overestimates  permeabilities.  Finally,  if  we  extrapolate  Shi¬ 
mizu’s  formula  to  our  case  we  calculate  a  permeability  of  2.5x  10"4  cm2 .  This  coincides  with 
Bader’s  upper  bound,  and  thus  is  taken  as  the  upper  limit  in  this  experiment.  A  lower  limit  is 
probably  in  the  range  of  1. Ox  1 0"4  cm2  since  Colbeck  and  Anderson’s  (1982)  data  suggest  a  per¬ 
meability  of  1 ,4x  10“4  cm2  for  our  snow  density. 

Fluid  properties  are  given  in  Table  4  and  a  thermal  conductivity  of  6.25x  1CT4  cal/cm  s  °C  is 
estimated  from  the  heat  flow  data.  Using  the  properties  of  air  at  -20°C,  the  lowest  critical  heat 
tlux  that  we  calculate  is  qcr  =  1 ,59x  10'3  cal/cm  s.  This  is  about  a  factor  of  6  greater  than  the  ob¬ 
served  onset  point,  a  large  discrepancy,  especially  for  a  lower  bound.  This  apparently  results  front 
a  large  error  in  the  calculated  permeability,  which  is  not  too  surprising  considering  the  uncertain¬ 
ty  in  permeability  discussed  above. 

The  temperature  profiles  shown  in  Figure  30  represent  two  cases  where  the  thermal  conductiv¬ 
ity  (Fig.  31 )  suggests  that  convection  is  occurring,  and  one  case  where  it  suggests  that  convection 
is  not  occurring.  The  lateral  profile  for  the  case  without  convection  in  Figure  30a  shows  a  drop 
near  the  boundaries,  but  this  may  simply  be  caused  by  a  slight  displacement  of  these  thermo¬ 
couples,  In  any  case  the  lateral  temperature  differences  are  too  small  compared  to  the  error  to 
draw  any  conclusions.  From  the  vertical  profiles  in  Figure  30b,  we  see  a  temperature  profile  that 
is  characteristic  of  the  centerline  of  a  single  longitudinal  roll  (see  Fig,  2).  One  roll  would  also  be 
expected  from  the  theoretical  results  of  Tewari  and  Torrance  (1981).  When  lateral  heat  fluxes  of 
the  order  observed  in  the  experiment  are  included  in  our  numerical  model,  we  find  that  only  very 
near  the  critical  Rayleigh  number  will  the  cellular  form  consist  of  two  units.  This  could  be  the 
case  for  the  one  experimental  run  shown  here  in  which  we  did  not  expect  convection.  This  Ray¬ 
leigh  number  is  the  one  nearest  the  estimated  onset  point  in  Figure  31 .  The  model  predicts  that 
close  to  the  critical  Rayleigh  number,  the  influence  of  convection  on  heat  transfer  is  minimal.  The 
presence  of  convection  is  hard  to  deduce  solely  front  the  temperature  profiles. 
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Theoretical  results  for  the  mode  of  convcc- 
—  tion  at  the  onset  are  not  available  for  the  con- 

-  0  65  /  v  /  /  stanl  ^ux  bottom  case.  However,  Tewari  and 

09 — &  40  T  7  f//  Torrance  ( 1981 )  solve  for  the  onset  form  of 

\  J  convection  with  an  isothermal  bottom  and 

vy  permeable  top.  From  Table  1  we  see  that  both 

^  06  —  the  critical  Rayleigh  number  and  the  critical 

9  M  wave  number  are  identical  for  the  cases  of  iso- 

y  thermal  bottom,  permeable  top  and  flux  bot- 

/  k  A  tom,  impermeable  top.  This  is  a  mathemat- 

/  y  J  J  ical  consequence  of  there  being  one  Neumann 

0  2  —  ^  yy //  type  (specific  gradient)  condition  for  each 

_  y  case.  We  *ben  expect  the  form  at  the  onset 

I  I  i  I  °f  convectl°n  t0  be  also  similar.  For  the  di- 

0  -0  2  -0  4  -0  6  -0  8  -io  mensions  of  the  current  experiment,  Tewari’s 

(t-t0)/at  results  indicated  that  four  hexagonal  cells 

should  be  present.  This  was  certainly  not  the 

Figure  28.  Vertical  temperature  profiles  for  case  with  the  present  experiment.  The  dis- 

convective  heat  transfer  in  glass  beads  and  agreement  may  be  ascribable  to  the  presence 

water.  Dashed  line  indicates  the  pure  con-  of  the  lateral  heat  fluxes  discussed  earlier. 

duction  solution.  Although  not  dominant  in  the  heat  transfer 

process,  they  may  be  significant  enough  to 
alter  the  form  of  convection. 

From  this  study  with  glass  beads,  three  conclusions  are  apparent.  First,  that  we  can  induce 
and  reasonably  measure  convection  in  porous  media.  Second,  that  it  is  necessary  to  better  con¬ 
trol  the  upper  thermal  boundary  and  the  heat  fluxes  through  the  sides  before  doing  experiments 
on  snow.  Third,  from  a  comparison  of  theoretical-numerical  results  with  the  experimental  re¬ 
sults,  we  can  reasonably  conclude  that  the  assumption  of  a  constant  flux  lower  boundary  condi¬ 
tion  was  a  sound  one. 


Figure  28.  Vertical  temperature  profiles  for 
convective  heat  transfer  in  glass  beads  and 
water.  Dashed  line  indicates  the  pure  con¬ 
duction  solution. 


Results  of  the  snow  experiment  are  shown  in  Figure  29.  Again,  effective  thermal  conductiv¬ 
ity  is  shown  as  a  function  of  heat  flux.  We  note  that  the  experimental  results  indicate  that  con- 
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Figure  29.  Effective  thermal  conduc 
tivitv  versus  net  vertical  heat  flux. 
Experimental  data  from  snow  ex¬ 
periments. 
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a.  y  direction  profile. 


b.  x  direction  profile. 


Figure  2  7.  Horizontal  temperature  profiles  for  convective  heat 
transfer  in  experiments  with  glass  beads  and  water. 


ing  down  the  sides.  The  magnitude  of  the  lateral  gradient  increases  as  the  Rayleigh  number  in¬ 
creases,  as  we  would  expect  for  such  a  cell. 

The  magnitude  of  the  lateral  gradients  could  be  enhanced  both  by  temperature  variations  on 
the  top  plate  and  by  heat  losses  out  the  sides.  The  center  of  the  top  plate  was  always  wanner 
than  the  edges,  and  thus  the  returning  fluid  lost  heat  as  it  passed  under  the  colder  outer  parts  of 
the  plate.  By  comparing  the  heat  fluxes  out  the  sides  with  the  heat  fluxes  that  would  have  oc¬ 
curred  if  the  lateral  profiles  were  only  attributable  to  conduction,  fluxes  out  the  sides  were  only 
about  10%  of  those  calculated  by  the  lateral  gradient,  and  thus  must  be  a  relatively  small  contrib¬ 
utor  to  the  lateral  temperature  profiles. 

The  three-dimensional  convection  observed  should  lead  to  enhanced  gradients  near  the  top 
and  bottom  surfaces,  and  to  nearly  isothermal  conditions  in  the  center.  Figure  28  shows  the  di¬ 
mensionless  centerline  temperature  profiles  for  the  three  convective  runs.  The  data  show  ihe  ex¬ 
pected  trend,  with  the  profile  becoming  more  distorted  as  convection  intensifies. 
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Table  5.  Properties  of  air  and  water. 


(dV/iT)p **  Ptt  P *** 

T(° C)  pjfg/cm3)*  p^g/cms)*  (S(° C'  )t  C C' )  (cal  s/cm1  (f  Cl  (cal  s/cms  ° C1 ) 


Air 


-20 

1.40X  10s 

1.61X  10* 

3.95X103 

- 

1.15X  10s 

_ 

-10 

1 .31  X  IO3 

1.67X  104 

3.80X103 

- 

0.98X  10s 

- 

0 

1.29X  103 

1.72X  104 

3.66X  103 

- 

0.85X  10s 

- 

10 

1.25X  10s 

1.76X  IO4 

3.53X  LO3 

- 

0.7SX105 

- 

20 

1.20X  103 

1.81X  104 

3.41  X  1 03 

- 

0.65X  10s 

_ 

40 

1.13X103 

1 .9 1 X  1 04 

3.19X  103 

- 

O.S1X  10s 

- 

60 

1.06X103 

2.00X  104 

3.00X1 03 

- 

0.40X  105 

- 

80 

1.00X  103 

2.09X  104 

3.83X  103 

- 

0.32X  10s 

- 

100 

0.95X  1 03 

2.17X  104 

2.68X  1 03 

- 

0.27X  10s 

Water 

10 

1.000X  103 

1.310X101 

_ 

10.37X  10! 

_ 

12.30X10* 

20 

0.998X  1 03 

1.000X  10! 

- 

2.00X  103 

_ 

1.99X10* 

30 

0.996X  1 03 

0.796X  103 

- 

2.9SX103 

— 

2.32X  10* 

40 

0.992X  103 

0.653X  103 

- 

3.80X  1 03 

- 

2.42X10* 

50 

0.988X  1 03 

0.547X103 

- 

4.57X103 

- 

2.41X  10* 

60 

0.983X  1 03 

0.466X  103 

- 

S.28X  103 

- 

2.34X10* 

70 

0.978X  103 

0.404  X  103 

~ 

5.96X103 

- 

2.25X  10* 

•Rcberson  and  Crowe  (1980) 


fHolman  (1976) 

••Dorsey  (1 940)-(3  V/bDp  =  0/p/ 


ff/*  = 


► *  *  p  - 


P 

,  dV 


get  qcr  =  6.2x  IO"2  cal/cm2  s  for  air  and  qcr  =  1 .3 lx  10"3  cal/cm2  s  for  water. 

Obviously,  qcr  for  air  is  well  out  of  the  range  of  the  experiments  and  thus  convection  is  not 
expected.  Since  the  range  of  physical  property  variations  with  temperature  is  small  compared  to 
the  difference  between  the  critical  and  experimental  heat  fluxes,  there  is  no  need  to  further  con¬ 
sider  the  temperature  dependence  of  air  properties.  The  discussion  of  the  glass  bead-air  experi¬ 
ment  is  thus  concluded,  having  shown  that  experiment  and  theory  agree  well. 

For  the  glass  bead  and  water  experiments,  choice  of  the  proper  temperature  is  important  be¬ 
cause  the  critical  heat  flux  for  the  onset  of  convection  predicted  theoretically  is  within  the  range 
observed  experimentally.  It  is  also  important  because  fluid  properties  vary  substantially,  espe¬ 
cially  at  low  temperatures.  The  runs  at  high  values  of  heat  flux  were  conducted  with  average  tem¬ 
peratures  in  the  25°  to  30°C  range,  outside  the  temperature  range  where  fluid  property  variations 
are  most  significant.  Total  temperature  differences  top  to  bottom  were  less  than  15°C.  Thus 
we  can  conclude  that  our  earlier  calculation  of  the  critical  heat  flux  is  within  the  range  of  5-10%. 
From  Figure  26b  we  see  that  our  calculated  critical  heat  flux  agrees  well  with  the  experimental 
results. 

Above  the  critical  Rayleigh  number  on  Figure  26b,  agreement  between  the  experiments  and 
the  prediction  of  the  numerical  model  are  fair.  Discrepancies  may  be  introduced  because  the 
model  is  two-dimensional,  while  either  two-dimensional  rolls  or  three-dimensional  hexagon- 
possible  in  the  experiments.  The  cell  form  induced  by  convection  is  indicated  by  the  steady- 
state  horizontal  temperature  profiles  at  the  vertical  midpoint.  Figure  27  shows  these  tempera¬ 
ture  profiles  for  the  data  points  in  which  convection  appeared  to  be  occurring.  Each  indicates 
that  a  three-dimensional  cell  existed,  with  warm  fluid  rising  in  the  center  and  cool  fluid  return- 
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a.  Data  from  experiments  are  compared  with  theoretical  predictions. 


b.  Theoretical-numerical  results  are  compared  with  experimental 
data  from  glass  beads-water  runs. 


Figure  26.  Effective  conductivity  versus  net  vertical  heat  flux  qnet. 


where 


K  = 


(pCJ 


p’f 


For  a  constant  flux  lower  boundary,  we  earlier  defined  A T as qH/km .  Thus  the  critical  heat 
flux  for  the  onset  of  convection  is  given  by 


Per  = 


P0(pCp)f@gKH 


(55) 


For  a  constant  flux  bottom  and  an  isothermal,  impermeable  top,  the  critical  Rayleigh  number 
is  (from  Table  1)  27.1  Using  km  and  K  from  Table  4  and  fluid  properties  from  Table  5,  we 


The  net  heat  flux  was  calculated  by  adding  all  heat  losses  or  gains  to  the  heat  input,  and  then 
dividing  by  the  area  available  for  vertical  heat  (low.  The  effective  thermal  conductivity  was  cal¬ 
culated  by  dividing  the  net  heat  flux  by  *he  average  gradient  across  the  entire  sample.  A  sample 
calculation  that  illustrates  this  proceH"  c  s  given  in  Appendix  C. 


EXPERIMENTAL  RESULTS  AND  DISCUSSION 
Glass  beads 

To  compare  the  numerical  results  with  the  experimental  results,  we  must  know  the  physical 
properties  of  each  fluid  and  porous  medium.  The  necessary  properties  of  the  bead  compact  are 
the  thermal  conductivity  km  and  the  intrinsic  permeability  K.  The  values  of  thermal  conductivi¬ 
ty  from  the  Luikov  et  al.  (1968)  model  are  in  close  agreement  with  those  measured  by  Bau  (1980), 
and  thus  that  model  is  a  useful  tool  in  our  work. 

Permeabilities  may  be  estimated  from  the  Carmen-Kozeny  relationship  (Wallis  1969) 


where  d  is  a  mean  particle  diameter  and  e  the  porosity.  Bau  (1980)  measured  permeabilities  for 
the  same  glass  beads  used  here  and  found  that  this  relation  describes  permeabilities  well.  Calcu¬ 
lated  parameters  for  our  experiments  are  listed  in  Table  4. 

Table  4  lists  temperature-dependent  fluid  properties  and  a  product  as  it  appears  in  the  Rayleigh 
number.  The  values  at  50°C  and  30°C  for  air  and  water  are  used  in  the  current  discussion  and, 
later,  the  effects  of  temperature-dependent  fluid  properties  are  discussed. 

The  essence  of  the  experimental  and  theoretical  results  with  glass  beads  is  shown  in  Figures  26a 
and  b,  which  are  plots  of  effective  thermal  conductivity  versus  net  heat  flux.  The  calculated  km 
corresponds  to  the  horizontal  line  in  Figure  26a,  and  in  both  cases  the  agreement  with  the  experi¬ 
ments  is  excellent.  The  inflection  point  in  Figure  26b  indicates  the  theoretical  onset  of  convection, 
with  a  corresponding  increase  in  the  rate  of  heat  transfer.  The  horizontal  lines  indicate  the  region 
where  the  heat  flux  is  less  than  the  critical  value  for  the  onset  of  convection.  The  value  of  ther¬ 
mal  conductivity  used  in  this  region  is  from  the  model  of  Luikov  et  al.  (1969)  and  is  tabulated  in 
Table  4.  Above  the  critical  heat  flux,  the  effective  thermal  conductivity  is  predicted  by  eq  44, 
which  is  a  result  from  our  numerical  model. 

The  critical  heat  flux  is  calculated  using  the  theory  of  convection  in  porous  media  presented  in 
the  Background-  Porous  Media  section.  The  Rayleigh  number  was  earlier  defined  as  (eq  5) 

pJgATHK 
Ra  = - 

I1K 


Table  4.  Physical  properties  for  glass  bead  experiments. 

Saturating  kf  km* 

fluid  (cal/cms^C)  e  (cal/cms°C)  Kf  (cm1  > 

Air  S.35X10'5  0.40  9.62XI0-4  2.53X10* 

Water  1,42X1 0'3  0.44  1,94X  10"*  3.77X10'* 

*  From  Luikov  et  al.  ( 1 968). 

fCalculated  from  Carmen-Kozeny  relation  (eq  54). 
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APPENDIX  A:  DERIVATION  OF  FINITE  DIFFERENCE  FORMULAE 
Constant  flux  bottom  (Fig.  Al)-no  phase  change 


:  -2 ; 

ifl  'i-1 _ ! _ :■*  I  j^l _ 

Figure  A1 .  Control  volume  at  lower  boundary. 
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Constant  flux  hot  tom -adiabatic  side  corner  (Fig.  A2) 


Figure  A  2.  Control  volume  at 
left  corner  of  lower  boundary. 

Accumulation  =  (T**/  ~ 


49 


n 


Conduction  =  - : - 

*x 

-  r  At  +  - 

Convection  =  -(t/7’)i+,/! 

h 

j  Ar  -  (wT)j 

Final  form: 

/npQ+i  'T'fi  \ 

V7  1  .1  _  7  1  ,1  > 

-  O 

A/ 

*2 

-  2 

Ax 

Opposite  corner,  final  form: 

lTi+1  _  rc  t 
t'nx.l  'nx,H 

A/  ‘  hi 

4.  2  +  a  (w7’)"* 

Formula  for  calculating  Nusselt  numbers  for 
isothermal  boundaries 

For  these  calculations,  we  assume  steady  state-bottom  boundary  (see  Fig.  Al): 


Accumulation  =  0 

.....  h 


Conduction 


K  A  .  (7U,,  -  r,tl)  h 

2  At  +  - K -  2 


K  A  *  (3T\  L 

2  A'  -  (Si,  *: 


h  h 

Convection  =  ( uT\_Vi  -  Ar  -  (uT)i+V2  —  At  -  (wT)l+)/t  hxAt. 


Final  form: 


(dT\  =  (7Yn,i  -  27M  +  ri-t,i)  K  +  (ri,2  ~  ri,i> 
laz/j  h*  2  +  hz 


hi  2 

K 


-  ((uT) j+vi  -  (M 2/j  “  Wmi  • 

Top  boundary,  final  form: 

an  (T’i+l.nz  -27i.nz  +ri-l.nz)  hz  (^i.nz-l  -^i.nz) 


Bottom  corners  (see  Fig.  A2): 


Accumulation  =  0 


Conduction  =  +  J'f)  h'At 

hx  2  /j2  2  \d2/,  2 


Convection  =  -(u7’)1  +  ,/j  yA/  -  (wD|t)4  jAr 


Final  form : 


(If)  =  ^.t  -  t'm)^  +  (t-i.2-  ^..t)  r  -  («n1+,,  -  wlt)i 


Opposite  corner,  Final  form: 


/0r\  hz  1 

fej  =  (T’nx-M  -  7-„x.i)  -  +  (rnx>2  -  rn8il)  r 


+  (ur)nx-W  Y  -  (w7Vv>- 


Top  corners.  Final  form,  i  =  1 : 


(92)  (7’2  nz  -  7",  nz)  -  (T,  ,nz-l  -  7*1  ,nz)  ^ 


+  («70i+h  7T  -  (*n. 


Top  corners,  final  form,  i  =  nx: 


l—\  =  _ (j  _  71  )  tj  _  j  ■»  J_ 

(02/  '  nx-l.nz  nx.nz-'  ^2  v  nx,nz-l  'nx.nz'  fj 


(m7\ x-V,  '  (^nz-A- 


The  Nusselt  number  is  then  given  by: 


Formula  for  calculating  Nusselt  number  for 
constant  flux  boundaries 


/  \ 

Nu  s  {(Tbot  -  7’top),v)-1  =(1  *.,(7-,.,  -riini))-‘  . 
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APPENDIX  B:  COMPUTER  PROGRAMS 


Main  program 


PROGRAM  70  CALCULATE  THERMAL  CONVECTION  PLOWS  IP!  FC*>CLS  PECIU» 
DIMENSION  ARRAYS 

CI*E NS  ION  T(75.21>.PSI(7E«21).SPAFE(75,Z1>,U( 75 «21)tYli Still 
Cl  PENSION  SOURCE  (75*21) 

Cl  PENS  I  ON  A(75)*fc(75),C(75)»R(75> 

CEFINE  ARRAYS 

t:  temperature 

PSl:  STREAM  FUNCTION 

SPARE:  SPARE  FIELD,  USEC  IN  CALCULATIONS  FOR  POTH  FSI  :NC  T 
u:  X  DIRECTION  VELOCITY  (HORIZONTAL) 

Ml  Y  DIRECTION  VELOCITY  (VERTICAL) 

SOURCE!  RIGHT  HANC  SICE  CF  MOMENTUM  EGUAT ION 
A.B.C  COEFFICIENTS  CF  IMPLICIT  FC  FCRM 
R  RIGHT  HAND  SIDE  OF  IMPLICIT  FD  FORM 

CEFINE  CONSTANTS 

CM:  SCR  RELAXATION  FACTOR 

HX:  GRID  SPACING  IN  X  DIRECTION 

hz:  GRID  SPACING  IN  Z  DIRECTION 

NX!  NUMBER  CF  NCCES  IN  X  TIRECTICN 
NZ:  NUMBER  OF  NODES  IN  Z  CIRECTICN 

ASRAT:  ASPECT  RATIO.  HOR  ./VERT . 

THETA:  ANGLE  CF  TILT,  FROM  MCRI2CNTAL 

ccnver:  ccnvergence  requirement  fcr  iterative 

SOLUTION  OF  STREAMFUNCT ICN  EQUATION. 

CELT:  TIME  STEP 

RA:  RAYLEIGH  NUMEER 

pert:  VALUE  CF  INITIAL  TEMPERATURE  DISTRIBUTION 
niter:  iteration  number  in  ADI  ITERATIVE 
NT:  TIME  STEP  NUMBER 

READ  CONSTANTS 

PRINT  70 
READ  (1,*)  RA 
PRINT  72 
READ  (i,»>  NX 
PRINT  82 
READ  <l,*l  NZ 
PRINT  83 

READ  <  1  •  *  >  ASRAT 
PRINT  8 A 

READ  (1,0  THETA 
PRINT  7A 
READ  (1,*)  CELT 
PRINT  75 

READ  (1,*)  IFILE 

IF  (IFILE. NE.l)  GOTO  15 

PRINT  76 

PRINT  77 

READ  <1,«>  FSEG 

PRINT  78 

RE  AO  (1,*)  IPRINT 

IF  (IPRINT  .NE.  1)  GOTO  16 

PRINT  81 

READ  (1,*)  PSE6 

PRINT  79 

READ  (1,*)  IPSI 

PRINT  83 

READ  (1.0  ITEMP 

FORMAT  ('INPUT-  SEQ:  (I.E.  IF  SEQ=5»EVERY  FIFTH  TIME  STFP') 

FORMAT  ('WILL  GENERATE  RESULTS  THAT  ARE  FILED)*) 

FORMAT  ('INPUT  IFRINT,  IPR1NT=1  VlLL  PRINT  RESULTS') 

FORMAT  ('INPUT  IFSI.  I  PS  I  =  1  GIVES  PERMEABLE  TCP') 

FORMAT  ('INPUT  ITEMP*  ITEMP=1  GIVES  CONSTANT  FLUX  POTTCM') 

FORMAT  ('INPUT  PSEQt  HOW  OFTEN  RESULTS  ARE  PRINTED') 

FORMAT  ('INPUT  NZ,  NUMBER  CF  NOOES  IN  Z  DIRECTION') 

FORMAT  ('INPUT  TFETA, ANGLE  OF  TILT  FROM  HOR.  IN  DEGREES') 

FORMAT  ('INPUT  ASPECT  RATIO  ( H OR  I ZO NT AL /VERT  I C AL ) • ) 

FORMAT  ('INPUT  OELTA  T»  TIME  STEP') 

FORMAT  ('INPUT  IFILE.  IF=1,  RESULTS  WILL  BE  FILED') 

FORMAT  ('INPUT  RAYLEIGH  NUMBER') 

FORMAT  ('NOTE  THAT  BOTH  X  AND  Z  NUMBERS  OF  NOOES') 

FORMAT  ('MUST  MATCH  THOSE  OF  THE  INPUT  FILES') 


rioooo  noo  oon 


FORMAT  (  ’INPUT  NX,  NUMBER  OF  NODE'  IN  X  DIRECTION*) 
FX=ASRAT/<NX-I> 

►  2  =  1  ./ <N2- 1 > 

AR=H2/HX 

THtTA=THETA.3.1N159/l&G. 

ARSQ=AR««2. 

N  T  =  L 

PER  T  =  1  .E-2 
N2M-NE-1 

N<M;N4-1 

C  M  =  1  •  7 

CONVER=l.f-A 
I  aRI  T  =  3 1 

CO  1C  v,  =  l,N2 

TLE  (u-l  > 

DO  1C  1=1, NX 
T  (  1  ,  u  )  =  T  L  E  V 
CONTINUE 

TEMP  PERTURBATION  TO  START  FLOWS 
T  (  NX*2)  =  T  (NX  ,2  MFERT 
START  CALCULATIONS  HERE 
MT  =  2 

IF  (HEMP. EG. 1)  NT  =  1 

M p  -  2M 

IF  (IPSI.EC.l)  MF=N2 

N  T  =N  T* 1 
PRINT  62, NT 
NI TER  =  C 

CZ=DELT/(2.»H2*«2> 

CX  =  C£LT/  (2  .«HX  •  *2  ) 

CALCULATE  SOURCE  TERM 

DO  20  U  =  2  »NZM 
CO  2  C  1=2.  NXN 

SCURCE(I,vi)iRA»AR»(T(I«l»vJ)-T(I-l,vJ>>»FZ»CCS(THETA)/2, 
r  n , TTnnr  -RA»P2*SIN(THETA>.{T<I,L»l)-Ttl,J-l>>/2. 

CUM  I  NUt 

SET  SPARE  FIELC  TC  2ER0  AT  BOUNDARIES 

FOR  SOLUTION  OF  S TR £ A MF U N C T I C N  EOLATION 

DO  21  J=1,NZ 

SP  A RE ( 1  * J >  =  C  •  C 
SPARE(NX»J)=0«0 
CONTINUE 
CO  22  1  =  1, NX 

SPARE (1,11=2.0 
SPARE (I »N  2  >  =  C  .  Q 
CONTINUE 

00  FOR  PERMEABLE  TOP  ONLY 

IF  (IPSI.NE.1I  GOTO  132 
DO  23  1=2, NXN 

SCURCE(l«NZ)=RA*(FZ/2.)*(CCS(TFETA)*»R*(T<l4l,NZ)-TU-l,NZ)> 

CONTINUE  *SIN(THETA)*(3.»T  (  I,N2)-A.*T (I,NZM)4T(I »  NZ  *  2  > )> 

AC  I  -SCR  STEFS  CALCULATE  PS  I  UNTIL  ADEQUATE  CONVERGENCE 

STEP  1,  IMPLICIT  IN  X 

MTER  =  MT£R*1 
OIFFMX=O.C 
PSIMAX=0.C 
PS  I  H1N  =  0 . 0 
CO  30  U=2,MP 

DO  32  1  =  2, NXM 
A ( I  -  1 ) =ARSC 
e<  I-l»=-2.*ARSQ-OM 
Cf I-1)=ARSC 
IF  (U.NE.N2)  GOTO  5 

R(I-1)=-2.«PSI(I,N2-1)4(2.-CM)*PSI(I,N2)*S0URCE tl,U> 

GOTO  32 

R(I-1)=SCURCE(I»J)-PSI(I,J»1)-PSI  (I.J-1)4(2.-CM)*PSI(I,J) 
CONTINUE 

CALL  THOMF  (A.6.C.R.NX-2) 

00  38  1=2  , NXM 

SPARE(I,u)=R(I“l) 

CONTINUE 

CONTINUE 

STEP  2,  IMFLICIT  IN  7 


34 

6 


35 

33 

C 

c 


120 

67 
C 

C  CALCULATE  VELOCITY  FIELO 

C 

00  36  I:2»NXN 

00  37  J:2,NZM 

U<I»«.)=-fPSI(I»J*l>-PSI<I.J-l>>/<2.*HZ> 

*<I»«.J  =  (FSI<l4i,JJ-PSI(I-l*i.>>/I2»»PX) 

37  CONTINUE 

U(I ,1 >=-(4. *PSI «I »2>-3.«PSI «I «1>  -PS  III »3) )/l2.*HZ> 
U<I,NZ>=-(-4.*PSl«l,NZHI*3.*PS]f  I»N2)*PSI(  I  ,  NZ- 2  »  >  /  (  2  .  *  F  2  ) 
36  CCMINLE 

CO  31  J=2,NZM 

y<l,J)  =  (4.«FSU2,J>-3.*PSIIl.J>-FSI(2,vJ>)/I2.*HX> 
b(NX,vi):(-4.«FSI(NXH,J)*3.*PSl<NX,J)+PSI(NX-2,U))/<;.*EX) 
31  CONTINUE 

C 

IF  IIPSI.NE.l)  GCTO  103 
CO  39  1:2, NXM 

U(I,NZ):<PSm*l,NZ)-PSl<  1-1,NE»)/(2.*HX) 

39  CONTINUE 

UIl,NZ):f4.*PSI<2»NZ>-3.*PSI ( 1 ,NZ > - FS II 3 ,NZ > > / < 2 . »H X > 

UINX,  NZJ:<3,*PSI(NX,N2)-4..PSHNXP,NZ)*PS  II  NX-2,  NZ  >>/«!.  *FX> 

C 

c  ncn  calculate  temperature  fielc 

C  USING  NON  ITERATIVE  ADI  AS  A  TIME  STEPPING  WEANS 

C 

C  IMPLICIT  IN  X 

C 

103  00  40  U:MT ,NZM 

DO  41  1:1,  NX 

IF  U.EQ.1J  GCTC  7 
WP=<W<I.J*l)*Uf I,J>>/2. 

WM:(ta<I,U-l)*y<I,v)))/2. 

R< I) :T< I ,U4l ).CZ»C2,-HZ*«WP-AESCUF) ») 

!  ♦T(I,fc>«(l«*CZ*(4,*PZ»(NF4A6S<UP>-VW*ABS(NP)>>) 

!  *T ( I  * U-l) *C2« <2.4HZ»IWH*ABS<WM  )> > 

GOTO  8 

7  WP:W«I,2)/2. 

R(I>:T<I,2»*CZ«2,*C2.-h2*ChF-ABStyP>>>  ♦ 

!  T (I ,1  )• I 1 ,-CZ* (4 .42  ,*H2*  (MP4ABSI WP) ) > >  * 

!  4,*CZ»FZ 

8  IF  (  I.EQ.l  )  GOTO  3 

IF  II.EO.HX)  GOTO  4 
EP:(UI«ln)4UIiJI/2. 

UH=fL(I-l,U)4U<I,J>>/2. 

A(I):CX4<-2.-HX4(UW4ABSIUM)l) 

e<I>:1.4CX«(4.4HX»(UP4APS<UF)-UW4AES<UW)») 

CII):CX»(-2.4pX»(UP-APSILP))> 

GOTO  41 

3  UP  =  U  (2  <  J  >  /  2  « 

Ea»:1.4CX.(4.42.»PX4(UF4AtS«LP))> 

C<1>=CX*<-4.«2,*HX*<0P-AES(IPM) 

GOTO  41 

4  UW  =  U  (NX-1  ,v)/2. 

A(NX)=CX*(-4.-2.4HX44UM4ABS(UM)n 

E(NX>:1.4CX4(4.-2.*HX*(LM-AES(UN)n 

41  CONTINUE 

95  FORMAT  <E12.4,3X,E12.4,3X,E12.4,3X,E12.4> 


CO  33  1  =  2, NXM 
00  34  J : 2 , w P 
oJ:U-l 
tUJ):l. 

e  uj>=-2  .-cm 

C<JJ)=1. 

R(Jd>:SCURCE(I,d)-ARSQ4<SPARE<l4l,vJ)4SPARE(I-l,J) 
-2  •  "SPARE  (  I  ,  U ) >-0M»SPARE  I  I  »  «. ) 

CONTINUE 
A (NZM  )  =  2. 

CALL  THCMF  I  A  ,e , C , R  , MP- 1 J 
OC  35  d:2,t*P 
JJ  =  d-l  . 

DIFF:ABSIR(JJ»-FSIU,J>» 

IF  (CIFF.CT.DIFFMXI  DIF  FMX : C  I  FF 
PS  I  (  I  ,«J>=R  ( v.U) 

IF  (FSI C  I  ,U> .GT  .PSINAX  )  PS  I RAX  =  PSI (I ,J> 

IF  IFSI C  I  ,u  )  .LT  .FSIMINJ  PS  I M I N  =  PS I  1 1 , U > 

CONTINUE 

CONTINUE 

CALCULATE  CONVERGENCE 

CRIT:0IFFMX/(AtS<PSlMAXj4AeS<FSIMIN)4;.'1} 

IF  ICR  IT  .G  1. CONNER)  GCTC  I '2 
PRINT  12C,  PSIMAX  ,PS  I  M  I  N 

FORMAT  (•PSIMAX:«,E12.4,3X,*PSIMIN:‘,E12.4) 

PRINT  61, NITER 

FORMAT  I ‘NITER:*, 14, *CRIT:*,E  12,4  > 


C 

c 

C 


SOLVE 


oon^n  non  ooo  non 


CALL  ThOMF  (A.B.C.R.NXJ 
00  43  1=1. NX 

SPARE ( 1  .  J  1  =  R  C  I  l 
96  FORMAT  (4X.E12.4) 

43  CONTINLE 

43  CONTINUE 

IMPLICIT  IN  2 

00  44  1  =  1  .NX 

SPARE (I*NZ)=*1.0 

44  CONTINUE 

00  53  1  =  1  . NX 

IF  (ITEMP.NE.il  GCTC  9 
VP=W(I,2>/2. 

B<l)  =  l.*C2«(4.*2.*HZ»<UP*A6SftiP>)> 
C(l>=C2«(-4..2.»hZ»(WP-ABS(WP)>> 

9  CC  51  *»  =  M  T  »  NZ  M 

JJ=J-MT*1 

IF  (U.EQ.ll  GOTO  11 
WP=( W( I,L«1 >«W( 1 ,U) I /2 . 

6M  =  (  k(  I, vl«W (  I.U-11  »  /2  . 

A(JJ)=CZ«(-2.-H2«(WM*ABS(VM>>> 

B  ( JJ  »=1.«CZ*I4.*HZ*<UP*ABS<  WP  l-VM*ABS( WM1  »  1 
CtJ«J)  =  C2*(-2.*HZ»CUP-ABSIUP>l  I 
11  IF  (I.EC.l)  GCTC  1 

IF  II.EQ.NX>  GOTO  2 
UP=(L(  I  .  J  >  «ll  (  I  *1  .  J  >  > /2  . 

UH= ( L ( I ,U)«Li I -1 , J I  I /2  . 

R«0L)  =  SPARMI»1.vJ>*CX.<2.-F>.  <UP-AES<UP>I>*SFARE  (  I  ,  .  I  « 
!  <1 .-CX» (4 ..HX» (UP*A8S (UP l-UM.ABS (UM >  1 1  I  ♦ 

!  SPARE (I -1 .J) »CX» 12 •  ♦  H  X  * (UM*ABS(UM)  >  ) 

GOTO  51 

1  UP=L(2«.)/2. 

R  < Jd)  =  SPARE(l»J>»(l.-CX«t4.42 . *HX « ( UP* ABS ( UP  I >>  >  * 

!  SPARE(2«JI«CX*(4.-2.»HX*(UP-AES(UP>)1 

GOTO  51 

2  LM  =  U  <N X -  1 . *  > /2 • 

R(JJ>=SPARE(NXM.J>»CX.(4.*2.*HX*(UM+ABS(UM>>>* 

!  SPARE(NX,JI«(l.*CX*(-4..2.*HX*(UM-ABS(yH>>>> 

51  CONTINUE 

F  CLLOW I NG  IS  FCR  CONSTANT  HEAT  FLUX  CNLt 

IF  (ITEMP.EQ.ll  R(ll=R(ll«4.«CZ*HZ 
R(NZ-MT)=R(N2-MTI-C(N2-HT)»SPARE(I.NZ> 

DO  53  J=MT  »NZM 
JJ=J-1 

IF  (MT.EC.ll  JJ=J 
53  CONTINUE 

SOLVE 

CALL  THCNF  ( A ,B «C , R »N Z -MT > 

DO  52  J=MT,N2M 

T(I.Ul=R(J-HT*l> 

52  CONTINUE 
50  CONTINUE 


NOU  PRINT  CR  FILE  IF  SO  OESIGNATEC 

IF  (IPRINT.NE.il  GOTO  I r 4 
AA  =  FLOAT (NT  )  /  FSEC 
bB=INT(NT/PSEC> 

IF  (A  A.  NE.B8  I  CCTC  1  "  4 
PRINT  65 

65  FORMAT  ( *  VELOCITY  FIELCM 

00  199  J  =  1  •  N  Z 
MiNZ*  1-U 

PRINT  B6.  W(2.M1 »W(4,MI  ,N(6»H  ,U(fc «M I ,V(  13  .Ml 

199  CONTINUE 
PRINT  66 

66  FORMAT  ( ‘TEMPERATURE  FIELD*) 

00  2u0  „=l»Ni 
M=NZ.1-J 

PRINT  88.  T(2»M)*T(4.M)»T(6»F1.T(P.P1,T(1'»M> 

200  CONTINUE 
C 

104  IF  (IF1LE.NE.1)  GCTO  135 
C 

AA=FL0AT(NT1/FSEC 
88  =  INT (NT/FSEQ  I 
IF  (AA.NE.EB)  GOTO  135 
PR  I NT  8  7 

87  FORMAT  (•  CO  YCL  WISH  TO  FILE,  YES=l.NO=C*> 

RE  AC  (l**l  IF  ILE2 
PRINT  187 

187  FORMAT  (*00  YOU  WISH  TO  CHANGE  TIME  STEP:  YES=1*1 
RC  AO  (1,*)  I T I  ME 
IF  (ITIME.NE.l)  CCTC  188 


nnnnnno 


PRINT  189 

189  FORMAT  ('INPUT  NEW  TIME  STEP'J 
RE  AC  lit*)  CELT 

188  PR  1  N  T  1 9  C 

190  FORMAT  ('DC  YOU  WISH  TO  CHANGE  RA,  YES=1'1 
RE  AO  (1,0  IRA 

IF  (IRA  .NE.  1)  CCTC  191 
PRINT  192 

192  FORMAT  ('INPUT  NEW  RA'l 
RE  AO  (1,0  R A 

191  PRINT  193 

193  FORMAT  OCC  YOU  WISP  TO  CHANCE  ANCLE.  YES=1*1 
READ  (1.0  IANG 

IF  (IANG.NE.il  GCTO  195 
PRINT  194 

1 9  A  FORMAT  ('INPUT  NEW  ANCLE.  IN  CEGRfESO 
READ  (1.0  ANG 
ANG=3.1A15?/18C. 

195  IF  (1FILE2.NE.1)  GOTO  1C5 
IWRITP  =  IWR  IT ♦ 1 

CALL  GFILEi  (IWRIT.2«'TFILENAME',C! 

CALL  GFILEJ  ( I W R I T P • 2  •  • PS  I F I L NM •  ,  B 1 
CO  lCjO  1=1. NX 

WRITE  (IWRIT.S91  (T(I.J)  .  J=1,N2) 

WRITE  (IWRITP.991  (PSKI.J)  •  ,.=  1,N2) 

1000  CONTINUE 

CALL  CLOSAL 
IWRIT=IWRIT.2 


NOW  COMPUTE  NUSSELT  NUMBERS  TOP  ANC  BOTTOM, 

ANC  CALCULATE  CCNVERGENCE  BASEC  CN  ENERGY 
CONSEHVATICN.  • 

BOTTOM  NUSSELT  NUMBER 
105  BOTNU=C.O 

00  59  1=2,  NXM 

UP=(U(I'I.1!«L(I, 111/2. 

UM  =  ( U  (  I  -1 ,1  ML  (  1  ,  11  1  /2. 

WP  =  W ( I , 2  I /2  . 

BOTNU  =  8CTNU-C.5*(AROT(I*1»1)*TII-1»1)“2,*T(  1,1)1)  - 
!  <T(  I  ,2  1  -  T  «  I  ,1  1  l/AR 

!  ♦ (HZ/A. )• ( T ( 1 ,1 )• (UP*ABS (UF)-UM*ABS(UM) 1 

!  *T(I*1,1)*(UP“APS(UF1)  -TC  I-l.l)*(UM*AfiS(UM)l) 

!  *(HX/2.)*((WP-A8S(WP))M  (I  ,2 1 ♦ (WP* ABS ( WP 1 1 »T  (  1,1)1 

59  CONTINUE 

UP  =  U  (2,11/2. 
wP=W(l,2)/2. 

BOTNU=BOTNL-C.E«(AR.(T<2,l)-T(l,ll).(T(l,2)-T(l.l))/AP) 

!  op* /a  .  i  •  <t  ( i ,  i  i*<up**esaiF ) } 

■  ♦  I  ( 2,  1 1  •  (LP-AeS( LP)  1  1  ♦  ( T (1  ,2).  IWP>ABS(WP) 1 

!  ♦  T(1,1)«(WP'*ABS(WP)))»(HX/A.  1 

UM=U(NX“l,l)/2. 

WP=w(NX,2)/2. 

OCTNU=BOTNL-O.E*(AR»(T(NXM,1)-T(NX,11)  ♦ 

!  (T (NX, 21-T (NX  ,  1 1  1 /AR  1  -  (HZ/A  .1  OT  (NXM  ,  1 1  * 

!  (UM*ABS(LM1)  *T (NX, 1 1  * (LM- A ES (UM)  1 1 

!  ♦(PX/A.)OT(NX,2)OWP-ABSUF>)*T(NX,:)» 

!  (WP*ABS( WP) 1 1 

c 

C  TOP  NUSSELT  NUMEER 

C 

T OPNU=  "  •  ’ 

00  fcj  1  =  2, NXM 

UP=(U(  I*1«N2)*U(I,N21) /2. 

UM=(UU-l,N2)*U(I,NZ))/2. 

W M  = ( W ( I  ,NZM),W(I,N2>)  /  2. 

TCPNU  =  TCPNU»C.:*AROT(I*l,N2),T(I-l,NZ)-2.*TU,NZ)l 
!  ♦(T(I,N2M)-T(I,N2))/AR 

!  - (HZ/ A. )• (T  < I ,N2) *( UP* ABS (UP »-UM*ABS(UN> 1 

!  *T(  I  ♦  1  ,NZ  1  »  (UP-ARS(UF  1 1  -  T(I-1,NZ1*(UMAABS(L'N))1 

!  ♦(PX/2.)*(T<1,NZ).(WM-AES(WM))*T(I,NZM)»(WM*AFS(WM)1) 

60  CONTINUE 
WM=(W(l,NZM)«W(l,NZ))/2. 

UP  =U ( 2 ,NZ 1/2. 

T0PNU=T0PNU,C.5OAR.(T(2,NZ)-T(l,NZ))*(T(l,NZM)-T(l,N2>)//F) 

!  -(HZ/A. ). (T(1,N2)«(UPaABS(UF  )  1 ♦ T ( 2 ,N 2 1 • ( UP - ABS (UF  1)1 

!  ♦(HX/A.).<T(1,NZM).(WM-ABS(WM))«T(1,NZ)»(WMAABS(WM))) 

WM=(W(NX,NZM)«W(NX,NZ))/2. 

UM=U(NXM,N2)/2. 

TOPNU=TOPNL*0 .5 • (AR» (T (NXM,NZ l-T (NX ,NZ)  1 1 
!  ♦0.5OT  (NX.NZM) -T(NX,NZ  )  l/AR 

•  *(P2/A.|.(T(NX,NZ)*(UM  ABS(LMl)  ♦  T(NXN,NZ1* 

!  (UM*  IBSCLM  1  1  1  ♦  (PX/ A  .  )  .  (T ( NX  ,NZ  )  * (WM- ABS  (  WM  1  1 

!  *T (NX.NZM) ,(WM«ABS(WM>  1  1 

TOPNU=TOPNU/ASRAT 
BOTNU=BOTNU/ASRAT 
BAL  =  (AES(UPNU-ECTNU)1/AES(EC1NU) 

PRINT  63,  TOPNU.EOTNU.BAL 
C 

IF  (NT.LT.100C)  GOTO  IOC 
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■jim."*.*ix.i,m.'X-"»-,l  mm" 


C 

C 

61  FORMAT  (IX, ‘NITERS*. If ) 

62  FORMAT  (f X,*NT=*,IC) 

63  FORMAT  (*TCPNU=*,1X,E12.A,5X,*FCTAU=*,1X,E12.A,IX, 

!  *PALANCE=*,1X,E12.A) 

88  FORMAT  <6 ( ZX  ,E  12  .A  )  ) 

99  FORMAT  <11  (€11. A, IX)) 

STOP 

END 

(INSERT  THOMFN 


Plotting  routine 
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C 

i  :o 

63 


75 

10 

80 

C 

c 

c 

c 
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65 

66 

67 

68 

69 

70 

C 

C 

c 

c 


PROGRAM  TC 
DIMENSION  ! 


PLOT  DATA  FILES 

(10C,21>,X(A>,2(A),rL<«),XPLT(‘>).ZPLT(4) 
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PRINT  50 

FORMAT  ('INPUT  TFETA, ANGLE  PET.EEN  X  AXIS  ANC  HORIZONTAL*) 
READ  (1.*)  THETA 
PRINT  51 

FORMAT  (‘INPUT  LENGTH  (IN  INCHES  )  OF  *.2  AXES*) 

RE  AC  (1,*)  XL.ZL 
CALL  PLOTS(C) 

ZMAX=1 . 

PRINT  52 
R  E  A  C  (  1 «  *  )  A  R 

FORMAT  (’INPUT  ASPECT  RATIO*) 

XM  A  X  =  A  R 
XLsXL/XMAX 
2L=Zl/2MAX 
PRINT  6  C 

FORMAT  (’INPUT  NUMBER  OF  HORIZONTAL  NODES*) 

R  E  A  0  <1,«)  MX 
FRINT  61 

FORMAT  ('INPUT  NUMBER  OF  VERTICAL  NCOES*) 

READ  (1,*)  M2 
C£LX=XMAX/ (MX- 1 ) 

CELZsZMAX/ (M2- 1 ) 

PRINT  63 

FORMAT  (‘INPUT  CATA  FILE  NAME*) 

CALL  GFILE1  (5,1, 'FILENAME*, 3) 

PRINT  75 

FORMAT  (  •  U  F  A  T  COLOR  PEN*) 

R  E  A  0  (1,*)  IPEN 
00  1C  1  =  1, MX 

READ  (5,80  (C(t,U),  J=  1 ,  M  2  1 
CONTINUE 
CALL  CLOSAL 
FORMAT  (11  (Ell. A, IX)) 

THIS  SECTION  SPECIFIES  THE  RANGE  ANO  CETAIL  CESIREC 
FOR  TFE  FLCT. 

DMAXsC  .0 
OMINsQ.O 
CO  2  C  1  =  1,  MX 
00  20  J  =  1  ,  **2 
IF  (0(1, J) 

IF  ( D ( I  ,  J  ) 

CONTINUE 
PRINT  65,  CM  A  X 

FORMAT  ( *  M  AX  V  A L L E  =  * , 1 X , E 1 2 . A > 

PRINT  66,  CHIN 

FORMAT  (  •  M  IN  V AL L E= ' , 1 X , E  1  2  .  A  ) 

PRINT  67 

FORMAT  (’INPUT  MAXIMUM  VALUE  CF  CCNTOUR  LINE  DESIRED*) 

READ  (1,»)  OMAX 
FRINT  68 

FORMAT  (’INPUT  CESIREO  MINIMUM  VALUE  CF  CONTOUR  LINES’) 
READ  (1,*)  OMIN 
PR i N T  69 

FORMAT  (’INPUT  CESIREC  NUMBER  CF  CONTOUR  LINES*) 

READ  (1,*)  NLINES 
PRINT  70 

FORMAT  (’INPUT  LINE  TYPE:  SEE  MANUAL  FOR  LINE  TYPE  CEF.’I 
RE  AO  (1,*)  ITYFE 

ACTUAL  PLOTTING  BEGINS  HERE 


.GT  .  DMAX  )  CM  A  X  =  0 (  I  ,  J  ) 
.LT.  CHIN)  CMINrOd.U) 


DELO=(CMAX-CMIN)/(NLltsES-l) 

DPLT=DH I N 

MXM=MX - 1 

MZM=MZ -1 

MsZMAX’ZL 
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CALL  N  £ UP  E  N (  IPEN) 

CO  25  N  =  l,  ALINES 
C  Q  *  3  l  =  1  .  M  *  M 
DO  30  J=1,M2M 

X(1)=CI-1)*DELX 
X ( 2  )  =  1  •CEL* 

X  <  3  >  -t  (  ?  ) 

x(4)=x( ;  ) 

1  <  1  )  =  (  1  >  *OE  L  2 

2 ( 3  >  =J*OEL2 

2  <  2  )  -  2  ( 1  ) 

2  (  4)  =  2  (  3  > 

CL(  1  >  =  D(  I  ,..) 

CL ( 2 )  =  D ( I O , J  ) 

CL ( 3  >  =  C  <  1*1* JO  ) 

CL(4  >  =  C (  I, v*l > 

I  P  L  T  - 1 
CO  31  K  =  l,4 
KF:K*1 

IF  (K*FC>4)  K  P  s  1 

IF  ((CPLT.GT.DL(K))  .AND  .  (DPLT.GT.DL(KP)))  G CTC  51 
IF  (ICPLT.LT.OL(K))  .AND.  <DPLT.LT.DL(KP>>)  C  TTC  31 
PCLATE=AeS((CPLT-DL(K>>/<CL(KP)-CL(K>>> 

XFL  T (  IPL  T  >  =  X (K  )  *XL 

IF  <2<K  ).F Q.2(KP>  )  XFLT (  1PLT  )  =  (X  (K  )  *POL ATE  * 

!  (  X  (  KF  > -X (K ) ) >  •  XL 

2  F  L  T  (  lf=LT>=Z€K»*ZL 

IF  (X(M).EQ.X(KP)  )  2FLT  f  1PLT >  =  (  2  (  K  )  ♦  P  0  L  A  T  E  * 

■  (2(KF)-2(K)J)*2L 

IPLT=IPLTO 
51  CONTINUE 

IF  (  IPLT  .EC  .  1  >  GOTO  2C 

IF  (IPLT.GT.3)  GOTO  30 

IF  (THETA.  EO.  01  GOTO  102 

CALL  ANGLE  (  X P L T (  1  )  , 2 P L T  ( 1 >  ,  T F E T A  ,  F  ) 

CALL  ANGIE  ( XPLT ( 2 > , 2PL T (2 >  ,T FE T  i  ,F ) 

102  CALL  PLOT  (XPLT ( 1  )  ,2PLT (  1) , *> 

CALL  PLOT  (XPLT  (2  )  ,2PLT  (2)  ,  IT YPE  ) 

30  CONTINUE 

CPLT=CFLT*CELC 
25  CONTINUE 
C 

PRINT  71 

71  FCRMAT  (*CC  YCL  U1SF  TC  PLOT  MCRE  DATA:  YES=1») 

READ  (1,0  NPLOT 

IF  (NPLOT  .EQ.  1)  GOTO  IK 
PRINT  72 

72  FORMAT  (  •  C  C  YCL  WISP  TO  PLCT  CRIC:  YES=1*> 

READ  (1,0  IGRID 

IF  (IGRID  .NE.  II  GOTO  121 
C 

C  PLOT  GRIC 

C 

PRINT  75 

READ  (1,0  IPEN 

CALL  NEUPEN(IPEN) 

DO  40  J=1,M2 
X1=0.0 

20=(J-1)*CEL2*2L 
22  =  20 

IF  (THETA. NE.C)  CALL  ANGLE  ( X 1  ,2 2  ,  THE T A , H  I 
CALL  PLCT  (XI, 22, 3) 

00  40  1  =2  ,  MX 

Xl=  (  1-1  »  .CELXOL 
22  =  2C 

IF  (THETA. NE.C)  CALL  ANCLE  <  X  1 , 22  »T  HE  T  A  ,H  > 

CALL  PLOT  (XI, 22, 2) 

40  CONTINUE 
C 

CO  41  1  =  1,  MX 

21  =  0. 

X 0  =  (  I  *  1  I  •  CEL  X  •  XL 
X  2=  XO 

IF  (THETA. NE.C)  CALL  ANGLE  (X 2  ,2 1 , THE T A , H  ) 

CALL  PLCT  ( *2,21,3) 

CC  41  «J  =  2, M2 

2l=(w*l)*CEL2»2L 
X2  =  XC 

IF  (THETA. NE.O)  CALL  A N C LE ( > 2  ,2  1  , T F E T A , F  ) 

CALL  PLOT  02,21,2) 

41  CONTINUE 
C 

C  PLOT  BOUNDARIES  IF  SO  OESIRED 
C 

101  PRINT  73 

73  FORMAT  («CC  TOE  USE  TO  PLCT  FCUNCARIFSI  YES=!*) 

READ  (1,*)  I BNC 

IF  (IBND.NE.l)  C-CTO  1'3 

PRINT  75 

READ  (1,0  IPEN 

CALL  NEMPf N ( IPEN) 
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DO  9 C  J=  1  *  2 

xuc.c 

IF  (J.EC.  1)  2 C=  ,.P 
IF  (J.EG.2)  ZC=ZL*Z*AX 
Z2  =  Z0 

r»i  i‘TcFJ?*?»?C^CHL  ANGL£  <X1  ,Z2,TFFTA,P> 
CALL  FLCT  (  X  i  ,  Z  c  ,  Z  ) 

DO  SC  1=2, HX 

Xl  =  <  I-1>*CELX*XL 

CALL^PL0T*(X1»Z2*2)*'  *NCLf  «  *  ’  •« ’^T  A  .1- > 
^92  CONTINUE 

CO  91  1=1,2 
Z1=0.C 

IF  (l.EG.l)  XC=C.C 
IF  < I .E  G . 2  )  XC  =  XL  *  XM AX 
Xi  =  XO 

IF  ITHtTA.NE.C)  CALL  ANGLE  <X2,Z1,THETA»H> 
CALL  PL  CT  (X2,Z1,2)  * 

DC  91  J=2,PZ 

Zl=tv>l)*CELZ«ZL 

Pall‘^[5t*^?;L^L  ANGLE  ,x2*Z1  , theta, H) 

91  CONTINUE 

C 

103  CALL  PLOT  <C.,:.,999> 

STOP 

END 

11NSERT  ANGLE 


APPENDIX  C:  SAMPLE  CALCULATIONS 


Data  from  snow  experiment  2  February  1984. 


Figure  Cl .  Relative  dimensions  of  snow  sample  along  with  measured 
temperatures  (°C). 


V  =  5.83  V 
R  =  47.4  n 

4>jn  =  I'2//?  =  0.71  W  =  0.1714  eal/s 
Side  fluxes  =  {I  ins)  (area  of  heat  transfer)  (A77Az) 
kins  is  the  thermal  u  d  ctivity  of  tfie  polystyrene  insulation,  and  is 
kins  =  6.2x  I0“s  eal/cm  s°C. 

Right  =  (6.2x  10  5 )  x  (50x  10  cm2  )x  (1.4°C/5cm)  =+8.47x10“-'  cal/s 

Left  =  (6.2x  10  s )  x  (50xl0cm2)x  (l.7°C/5cm)  =  +1 ,03x  10'2  cal/s 

F  ront  =  (6.2x  I0“5 1  x  (50x40  cm2 1  x  (l  ,0°C/10  cm)  =  -1 .21  x  10~2  cal/s 

Back  =  (6.2x  10“s )  x  (50x40cm2)x  (1 .4°C/I0 cm)  =  -1  69x  10'2  cal/s 

sum  of  side  fluxes  =  -1.02x  10~2  cal/s 

",  of  input  =  0.0102/0.1714  =  5 .99? 

Bottom  =  (6.2x  10s )  x  (40x  10cm2)x  (5.1°C/1 0  cm)=  -1 .26x  10"2  cal/s 

Net  heat  flow  =  0.1 714 

-0.0102 
-0.0126 

0.1486  cal/s 

net  heat  flux  =  net  heat  flow  per  area 

=  0.1486  cal/s  per  400  cm2 
<7ne(  =  3.7x  10“4cal/cm2s 

effective  thermal  conductivity  (fceff)  =  qnetH&TI&z) 

Aofr  =  3.7x  1 0-4 /((  —  4.8— ( -20. 1  ))/50 ) 

*eff  =  1 .18x  I0~'  cal/cm  s°C. 

Uncertainty  in  keff  is  estimated  by  multiplying  fraction  heat  loss  through  the  sides  by  the  ther¬ 
mal  conductivity. 

uncertainty  in  *eff  =  (1.18x  10’’)  x  0.059  =  ±6.96x  10‘5  cal/cm  s°C. 


A  facsimile  catalog  card  in  Library  of  Congress  MARC 
format  is  reproduced  below. 
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